Actual source code: ex27.c
1: static char help[] = "Poisson Problem in 2d and 3d with simplicial finite elements in both primal and mixed form.\n\
2: We solve the Poisson problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
4: This example solves mixed form equation to get the flux field to calculate flux norm. We then use that for adaptive mesh refinement. \n\n\n";
6: /*
7: The primal (or original) Poisson problem, in the strong form, is given by,
9: \begin{align}
10: - \nabla \cdot ( \nabla u ) = f
11: \end{align}
12: where $u$ is potential.
14: The weak form of this equation is
16: \begin{align}
17: < \nabla v, \nabla u > - < v, \nabla u \cdot \hat{n} >_\Gamma - < v, f > = 0
18: \end{align}
20: The mixed Poisson problem, in the strong form, is given by,
22: \begin{align}
23: q - \nabla u &= 0 \\
24: - \nabla \cdot q &= f
25: \end{align}
26: where $u$ is the potential and $q$ is the flux.
28: The weak form of this equation is
30: \begin{align}
31: < t, q > + < \nabla \cdot t, u > - < t \cdot \hat{n}, u >_\Gamma &= 0 \\
32: <v, \nabla \cdot q> - < v, f > &= 0
33: \end{align}
35: We solve both primal and mixed problem and calculate the error in the flux norm, namely || e || = || q^m_h - \nabla u^p_h ||. Here superscript 'm' represents field from mixed form and 'p' represents field from the primal form.
37: The following boundary conditions are prescribed.
39: Primal problem:
40: \begin{align}
41: u = u_0 on \Gamma_D
42: \nabla u \cdot \hat{n} = g on \Gamma_N
43: \end{align}
45: Mixed problem:
46: \begin{align}
47: u = u_0 on \Gamma_D
48: q \cdot \hat{n} = g on \Gamma_N
49: \end{align}
50: __________\Gamma_D_____________
51: | |
52: | |
53: | |
54: \Gamma_N \Gamma_N
55: | |
56: | |
57: | |
58: |_________\Gamma_D_____________|
60: To visualize the automated adaptation
62: -dm_adapt_pre_view draw -dm_adapt_view draw -draw_pause -1 -geometry 0,0,1024,1024
64: and to compare with a naice gradient estimator use
66: -adaptor_type gradient
68: To see a sequence of adaptations
70: -snes_adapt_sequence 8 -adaptor_monitor_error draw::draw_lg
71: -dm_adapt_pre_view draw -dm_adapt_iter_view draw -dm_adapt_view draw -draw_pause 1 -geometry 0,0,1024,1024
73: To get a better view of the by-hand process, use
75: -dm_view hdf5:${PWD}/mesh.h5
76: -primal_sol_vec_view hdf5:${PWD}/mesh.h5::append
77: -flux_error_vec_view hdf5:${PWD}/mesh.h5::append
78: -exact_error_vec_view hdf5:${PWD}/mesh.h5::append
79: -refine_vec_view hdf5:${PWD}/mesh.h5::append
80: -adapt_dm_view draw -draw_pause -1
82: This is also possible with the automated path
84: -dm_view hdf5:${PWD}/mesh.h5
85: -adapt_primal_sol_vec_view hdf5:${PWD}/mesh.h5::append
86: -adapt_error_vec_view hdf5:${PWD}/mesh.h5::append
87: -adapt_vec_view hdf5:${PWD}/mesh.h5::append
88: */
90: #include <petsc/private/petscfeimpl.h>
91: #include <petscdmplex.h>
92: #include <petscsnes.h>
93: #include <petscdmadaptor.h>
94: #include <petscds.h>
95: #include <petscviewerhdf5.h>
96: #include <petscbag.h>
98: PETSC_EXTERN PetscErrorCode SetupMixed(DMAdaptor, DM);
100: typedef enum {
101: SOL_QUADRATIC,
102: SOL_TRIG,
103: SOL_SENSOR,
104: SOL_UNKNOWN,
105: NUM_SOL_TYPE
106: } SolType;
107: const char *SolTypeNames[NUM_SOL_TYPE + 4] = {"quadratic", "trig", "sensor", "unknown", "SolType", "SOL_", NULL};
109: typedef struct {
110: PetscBag param;
111: SolType solType;
112: PetscBool byHand;
113: PetscInt numAdapt;
114: } AppCtx;
116: typedef struct {
117: PetscReal k;
118: } Parameter;
120: /* Exact solution: u = x^2 + y^2 */
121: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
122: {
123: PetscInt d;
125: u[0] = 0.0;
126: for (d = 0; d < dim; ++d) u[0] += x[d] * x[d];
127: return PETSC_SUCCESS;
128: }
129: /* Exact solution: q = (2x, 2y) */
130: static PetscErrorCode quadratic_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
131: {
132: PetscInt c;
133: for (c = 0; c < Nc; ++c) u[c] = 2.0 * x[c];
134: return PETSC_SUCCESS;
135: }
137: /* Exact solution: u = sin( n \pi x ) * sin( n \pi y ) */
138: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
139: {
140: const PetscReal n = 5.5;
142: u[0] = 1.0;
143: for (PetscInt d = 0; d < dim; ++d) u[0] *= PetscSinReal(n * PETSC_PI * x[d]);
144: return PETSC_SUCCESS;
145: }
146: static PetscErrorCode trig_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
147: {
148: const PetscReal n = 5.5;
150: for (PetscInt c = 0; c < Nc; c++) u[c] = n * PETSC_PI * PetscCosReal(n * PETSC_PI * x[c]) * PetscSinReal(n * PETSC_PI * x[Nc - c - 1]);
151: return PETSC_SUCCESS;
152: }
154: /*
155: Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation:
157: f:[-1, 1]x[-1, 1] \to R,
158: f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy)
160: (mapped to have domain [0,1] x [0,1] in this case).
161: */
162: static PetscErrorCode sensor_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
163: {
164: const PetscReal xref = 2. * x[0] - 1.;
165: const PetscReal yref = 2. * x[1] - 1.;
166: const PetscReal xy = xref * yref;
168: u[0] = PetscSinReal(50. * xy);
169: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) u[0] *= 0.01;
171: return PETSC_SUCCESS;
172: }
174: /* Flux is (cos(50xy) * 50y/100, cos(50xy) * 50x/100) if |xy| > 2\pi/50 else (cos(50xy) * 50y, cos(50xy) * 50x) */
175: static PetscErrorCode sensor_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
176: {
177: const PetscReal xref = 2. * x[0] - 1.;
178: const PetscReal yref = 2. * x[1] - 1.;
179: const PetscReal xy = xref * yref;
181: u[0] = 50. * yref * PetscCosReal(50. * xy) * 2.0;
182: u[1] = 50. * xref * PetscCosReal(50. * xy) * 2.0;
183: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) {
184: u[0] *= 0.01;
185: u[1] *= 0.01;
186: }
187: return PETSC_SUCCESS;
188: }
190: /* We set up residuals and Jacobians for the primal problem. */
191: static void f0_quadratic_primal(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[])
192: {
193: f0[0] = 4.0;
194: }
196: static void f0_trig_primal(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[])
197: {
198: const PetscReal n = 5.5;
200: f0[0] = -2.0 * PetscSqr(n * PETSC_PI) * PetscSinReal(n * PETSC_PI * x[0]) * PetscSinReal(n * PETSC_PI * x[1]);
201: }
203: static void f0_sensor_primal(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[])
204: {
205: const PetscReal xref = 2. * x[0] - 1.;
206: const PetscReal yref = 2. * x[1] - 1.;
207: const PetscReal xy = xref * yref;
209: f0[0] = -2500.0 * PetscSinReal(50. * xy) * (xref * xref + yref * yref) * 4.0;
210: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) f0[0] *= 0.01;
211: }
213: static void f1_primal(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[])
214: {
215: const PetscReal k = PetscRealPart(constants[0]);
217: for (PetscInt d = 0; d < dim; ++d) f1[d] = k * u_x[uOff_x[0] + d];
218: }
220: static void f0_quadratic_bd_primal(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
221: {
222: const PetscReal k = PetscRealPart(constants[0]);
223: PetscScalar flux;
225: PetscCallAbort(PETSC_COMM_SELF, quadratic_q(dim, t, x, dim, &flux, NULL));
226: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux * n[d];
227: }
229: static void f0_trig_bd_primal(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
230: {
231: const PetscReal k = PetscRealPart(constants[0]);
232: PetscScalar flux;
234: PetscCallAbort(PETSC_COMM_SELF, trig_q(dim, t, x, dim, &flux, NULL));
235: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux * n[d];
236: }
238: static void f0_sensor_bd_primal(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
239: {
240: const PetscReal k = PetscRealPart(constants[0]);
241: PetscScalar flux[2];
243: PetscCallAbort(PETSC_COMM_SELF, sensor_q(dim, t, x, dim, flux, NULL));
244: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux[d] * n[d];
245: }
247: static void g3_primal_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[])
248: {
249: const PetscReal k = PetscRealPart(constants[0]);
251: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = k;
252: }
254: /* Now we set up the residuals and Jacobians mixed problem. */
255: static void f0_mixed_quadratic_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[])
256: {
257: f0[0] = 4.0;
258: for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d];
259: }
260: static void f0_mixed_trig_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[])
261: {
262: PetscReal n = 5.5;
264: f0[0] = -2.0 * PetscSqr(n * PETSC_PI) * PetscSinReal(n * PETSC_PI * x[0]) * PetscSinReal(n * PETSC_PI * x[1]);
265: for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d];
266: }
267: static void f0_mixed_sensor_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[])
268: {
269: const PetscReal xref = 2. * x[0] - 1.;
270: const PetscReal yref = 2. * x[1] - 1.;
271: const PetscReal xy = xref * yref;
273: f0[0] = -2500.0 * PetscSinReal(50. * xy) * (xref * xref + yref * yref) * 4.0;
274: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) f0[0] *= 0.01;
275: for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d];
276: }
278: static void f0_mixed_q(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[])
279: {
280: for (PetscInt d = 0; d < dim; d++) f0[d] = u[uOff[0] + d];
281: }
283: static void f1_mixed_q(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[])
284: {
285: const PetscReal k = PetscRealPart(constants[0]);
287: for (PetscInt d = 0; d < dim; d++) f1[d * dim + d] = k * u[uOff[1]];
288: }
290: static void f0_quadratic_bd_mixed_q(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
291: {
292: const PetscReal k = PetscRealPart(constants[0]);
293: PetscScalar potential;
295: PetscCallAbort(PETSC_COMM_SELF, quadratic_u(dim, t, x, dim, &potential, NULL));
296: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * potential * n[d];
297: }
299: static void f0_trig_bd_mixed_q(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
300: {
301: const PetscReal k = PetscRealPart(constants[0]);
302: PetscScalar potential;
304: PetscCallAbort(PETSC_COMM_SELF, trig_u(dim, t, x, dim, &potential, NULL));
305: for (PetscInt d = 0; d < dim; ++d) f0[d * dim + d] = -k * potential * n[d];
306: }
308: static void f0_sensor_bd_mixed_q(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
309: {
310: const PetscReal k = PetscRealPart(constants[0]);
311: PetscScalar potential;
313: PetscCallAbort(PETSC_COMM_SELF, sensor_u(dim, t, x, dim, &potential, NULL));
314: for (PetscInt d = 0; d < dim; ++d) f0[d * dim + d] = -k * potential * n[d];
315: }
317: /* <v, \nabla\cdot q> */
318: static void g1_mixed_uq(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 g1[])
319: {
320: for (PetscInt d = 0; d < dim; d++) g1[d * dim + d] = -1.0;
321: }
323: /* < t, q> */
324: static void g0_mixed_qq(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 g0[])
325: {
326: for (PetscInt d = 0; d < dim; d++) g0[d * dim + d] = 1.0;
327: }
329: /* <\nabla\cdot t, u> = <\nabla t, Iu> */
330: static void g2_mixed_qu(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 g2[])
331: {
332: const PetscReal k = PetscRealPart(constants[0]);
334: for (PetscInt d = 0; d < dim; d++) g2[d * dim + d] = k;
335: }
337: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *user)
338: {
339: PetscFunctionBeginUser;
340: PetscOptionsBegin(comm, "", "Flux norm error in primal poisson problem Options", "DMPLEX");
341: user->byHand = PETSC_TRUE;
342: user->numAdapt = 1;
343: user->solType = SOL_QUADRATIC;
345: PetscCall(PetscOptionsGetBool(NULL, NULL, "-by_hand", &user->byHand, NULL));
346: PetscCall(PetscOptionsGetInt(NULL, NULL, "-num_adapt", &user->numAdapt, NULL));
347: PetscCall(PetscOptionsEnum("-sol_type", "Type of exact solution", "ex27.c", SolTypeNames, (PetscEnum)user->solType, (PetscEnum *)&user->solType, NULL));
348: PetscOptionsEnd();
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: static PetscErrorCode SetupParameters(PetscBag bag, AppCtx *user)
353: {
354: Parameter *param;
356: PetscFunctionBeginUser;
357: PetscCall(PetscBagGetData(bag, (void **)¶m));
358: PetscCall(PetscBagSetName(bag, "par", "Poisson parameters"));
359: PetscCall(PetscBagRegisterReal(bag, ¶m->k, 1.0, "k", "Thermal conductivity"));
360: PetscCall(PetscBagSetFromOptions(bag));
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
365: {
366: PetscFunctionBeginUser;
367: PetscCall(DMCreate(comm, dm));
368: PetscCall(DMSetType(*dm, DMPLEX));
369: PetscCall(DMSetFromOptions(*dm));
370: PetscCall(DMSetApplicationContext(*dm, &user));
371: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
376: {
377: PetscDS ds;
378: DMLabel label;
379: PetscInt id, bd;
380: Parameter *param;
381: PetscWeakForm wf;
383: PetscFunctionBeginUser;
384: PetscCall(DMGetDS(dm, &ds));
385: PetscCall(DMGetLabel(dm, "marker", &label));
387: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_primal_uu));
389: switch (user->solType) {
390: case SOL_QUADRATIC:
391: PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_primal, f1_primal));
393: id = 1;
394: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))quadratic_u, NULL, user, NULL));
396: id = 2;
397: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
398: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
399: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_primal, 0, NULL));
401: id = 3;
402: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))quadratic_u, NULL, user, NULL));
404: id = 4;
405: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
406: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
407: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_primal, 0, NULL));
409: PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_u, user));
410: break;
411: case SOL_TRIG:
412: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_primal, f1_primal));
414: id = 1;
415: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL));
417: id = 2;
418: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
419: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
420: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_primal, 0, NULL));
422: id = 3;
423: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL));
425: id = 4;
426: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
427: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
428: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_primal, 0, NULL));
430: PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user));
431: break;
432: case SOL_SENSOR:
433: PetscCall(PetscDSSetResidual(ds, 0, f0_sensor_primal, f1_primal));
435: id = 1;
436: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))sensor_u, NULL, user, NULL));
438: id = 2;
439: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
440: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
441: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_primal, 0, NULL));
443: id = 3;
444: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))sensor_u, NULL, user, NULL));
446: id = 4;
447: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
448: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
449: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_primal, 0, NULL));
451: PetscCall(PetscDSSetExactSolution(ds, 0, sensor_u, user));
452: break;
453: default:
454: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]);
455: }
457: /* Setup constants */
458: {
459: PetscCall(PetscBagGetData(user->param, (void **)¶m));
460: PetscScalar constants[1];
462: constants[0] = param->k;
464: PetscCall(PetscDSSetConstants(ds, 1, constants));
465: }
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: static PetscErrorCode SetupPrimalDiscretization(DM dm, AppCtx *user)
470: {
471: DM cdm = dm;
472: PetscFE fe[1];
473: DMPolytopeType ct;
474: PetscInt dim, cStart;
475: MPI_Comm comm;
477: PetscFunctionBeginUser;
478: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
479: PetscCall(DMGetDimension(dm, &dim));
481: /* Create finite element */
482: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
483: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
484: PetscCall(PetscFECreateByCell(comm, dim, 1, ct, "primal_potential_", PETSC_DEFAULT, &fe[0]));
485: PetscCall(PetscObjectSetName((PetscObject)fe[0], "primal potential"));
487: /* Set discretization and boundary conditions for each mesh */
488: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
489: PetscCall(DMCreateDS(dm));
490: while (cdm) {
491: PetscCall(DMCopyDisc(dm, cdm));
492: PetscCall(DMGetCoarseDM(cdm, &cdm));
493: }
495: PetscCall(PetscFEDestroy(&fe[0]));
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: static PetscErrorCode SetupMixedProblem(DM dm, AppCtx *user)
500: {
501: PetscDS ds;
502: DMLabel label;
503: PetscInt id, bd;
504: Parameter *param;
505: PetscWeakForm wf;
507: PetscFunctionBeginUser;
508: PetscCall(DMGetDS(dm, &ds));
509: PetscCall(DMGetLabel(dm, "marker", &label));
511: /* Residual terms */
512: PetscCall(PetscDSSetResidual(ds, 0, f0_mixed_q, f1_mixed_q));
514: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_mixed_qq, NULL, NULL, NULL));
515: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_mixed_qu, NULL));
517: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_mixed_uq, NULL, NULL));
519: switch (user->solType) {
520: case SOL_QUADRATIC:
521: PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_quadratic_u, NULL));
523: id = 1;
524: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
525: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
526: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_mixed_q, 0, NULL));
528: id = 2;
529: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 0, 0, NULL, (void (*)(void))quadratic_q, NULL, user, NULL));
531: id = 4;
532: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 0, 0, NULL, (void (*)(void))quadratic_q, NULL, user, NULL));
534: id = 3;
535: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
536: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
537: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_mixed_q, 0, NULL));
539: PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_q, user));
540: PetscCall(PetscDSSetExactSolution(ds, 1, quadratic_u, user));
541: break;
542: case SOL_TRIG:
543: PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_trig_u, NULL));
545: id = 1;
546: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
547: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
548: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_mixed_q, 0, NULL));
550: id = 2;
551: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (void (*)(void))trig_q, NULL, user, NULL));
553: id = 4;
554: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (void (*)(void))trig_q, NULL, user, NULL));
556: id = 3;
557: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
558: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
559: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_mixed_q, 0, NULL));
561: PetscCall(PetscDSSetExactSolution(ds, 0, trig_q, user));
562: PetscCall(PetscDSSetExactSolution(ds, 1, trig_u, user));
563: break;
564: case SOL_SENSOR:
565: PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_sensor_u, NULL));
567: id = 1;
568: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
569: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
570: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_mixed_q, 0, NULL));
572: id = 2;
573: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (void (*)(void))sensor_q, NULL, user, NULL));
575: id = 4;
576: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (void (*)(void))sensor_q, NULL, user, NULL));
578: id = 3;
579: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
580: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
581: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_mixed_q, 0, NULL));
583: PetscCall(PetscDSSetExactSolution(ds, 0, sensor_q, user));
584: PetscCall(PetscDSSetExactSolution(ds, 1, sensor_u, user));
585: break;
586: default:
587: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]);
588: }
590: /* Setup constants */
591: {
592: PetscCall(PetscBagGetData(user->param, (void **)¶m));
593: PetscScalar constants[1];
595: constants[0] = param->k;
597: PetscCall(PetscDSSetConstants(ds, 1, constants));
598: }
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: static PetscErrorCode SetupMixedDiscretization(DM dm, AppCtx *user)
603: {
604: DM cdm = dm;
605: PetscFE fe[2];
606: DMPolytopeType ct;
607: PetscInt dim, cStart;
608: MPI_Comm comm;
610: PetscFunctionBeginUser;
611: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
612: PetscCall(DMGetDimension(dm, &dim));
614: /* Create finite element */
615: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
616: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
617: PetscCall(PetscFECreateByCell(comm, dim, dim, ct, "mixed_flux_", PETSC_DEFAULT, &fe[0]));
618: PetscCall(PetscObjectSetName((PetscObject)fe[0], "mixed flux"));
619: /* NOTE: Set the same quadrature order as the primal problem here or use the command line option. */
621: PetscCall(PetscFECreateByCell(comm, dim, 1, ct, "mixed_potential_", PETSC_DEFAULT, &fe[1]));
622: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
623: PetscCall(PetscObjectSetName((PetscObject)fe[1], "mixed potential"));
625: /* Set discretization and boundary conditions for each mesh */
626: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
627: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
628: PetscCall(DMCreateDS(dm));
629: while (cdm) {
630: PetscCall(DMCopyDisc(dm, cdm));
631: PetscCall(DMGetCoarseDM(cdm, &cdm));
632: }
633: PetscCall(PetscFEDestroy(&fe[0]));
634: PetscCall(PetscFEDestroy(&fe[1]));
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: PetscErrorCode SetupMixed(DMAdaptor adaptor, DM mdm)
639: {
640: AppCtx *ctx;
642: PetscFunctionBeginUser;
643: PetscCall(DMGetApplicationContext(mdm, (void **)&ctx));
644: PetscCall(SetupMixedDiscretization(mdm, ctx));
645: PetscCall(SetupMixedProblem(mdm, ctx));
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: int main(int argc, char **argv)
650: {
651: DM dm, mdm; /* problem specification */
652: SNES snes, msnes; /* nonlinear solvers */
653: Vec u, mu; /* solution vectors */
654: Vec fluxError, exactError; /* Element wise error vector */
655: PetscReal fluxNorm, exactNorm; /* Flux error norm */
656: AppCtx user; /* user-defined work context */
658: PetscFunctionBeginUser;
659: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
660: PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.param));
661: PetscCall(SetupParameters(user.param, &user));
662: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
664: // Set up and solve primal problem
665: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
666: PetscCall(DMSetApplicationContext(dm, &user));
667: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
668: PetscCall(SNESSetDM(snes, dm));
670: PetscCall(SetupPrimalDiscretization(dm, &user));
671: PetscCall(SetupPrimalProblem(dm, &user));
672: PetscCall(DMCreateGlobalVector(dm, &u));
673: PetscCall(VecSet(u, 0.0));
674: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
675: PetscCall(SNESSetFromOptions(snes));
676: PetscCall(DMSNESCheckFromOptions(snes, u));
678: for (PetscInt a = 0; a < user.numAdapt; ++a) {
679: if (a > 0) {
680: char prefix[16];
682: PetscCall(PetscSNPrintf(prefix, 16, "a%d_", (int)a));
683: PetscCall(SNESSetOptionsPrefix(snes, prefix));
684: }
685: PetscCall(SNESSolve(snes, NULL, u));
687: // Needed if you allow SNES to refine
688: PetscCall(SNESGetSolution(snes, &u));
689: PetscCall(VecGetDM(u, &dm));
690: }
692: PetscCall(PetscObjectSetName((PetscObject)u, "Primal Solution "));
693: PetscCall(VecViewFromOptions(u, NULL, "-primal_sol_vec_view"));
695: if (user.byHand) {
696: // Set up and solve mixed problem
697: PetscCall(DMClone(dm, &mdm));
698: PetscCall(SNESCreate(PETSC_COMM_WORLD, &msnes));
699: PetscCall(SNESSetOptionsPrefix(msnes, "mixed_"));
700: PetscCall(SNESSetDM(msnes, mdm));
702: PetscCall(SetupMixedDiscretization(mdm, &user));
703: PetscCall(SetupMixedProblem(mdm, &user));
704: PetscCall(DMCreateGlobalVector(mdm, &mu));
705: PetscCall(VecSet(mu, 0.0));
706: PetscCall(DMPlexSetSNESLocalFEM(mdm, PETSC_FALSE, &user));
707: PetscCall(SNESSetFromOptions(msnes));
709: PetscCall(DMSNESCheckFromOptions(msnes, mu));
710: PetscCall(SNESSolve(msnes, NULL, mu));
711: PetscCall(PetscObjectSetName((PetscObject)mu, "Mixed Solution "));
712: PetscCall(VecViewFromOptions(mu, NULL, "-mixed_sol_vec_view"));
714: // Create the error space of piecewise constants
715: DM dmErr;
716: PetscFE feErr;
717: DMPolytopeType ct;
718: PetscInt dim, cStart;
720: PetscCall(DMClone(dm, &dmErr));
721: PetscCall(DMGetDimension(dmErr, &dim));
722: PetscCall(DMPlexGetHeightStratum(dmErr, 0, &cStart, NULL));
723: PetscCall(DMPlexGetCellType(dmErr, cStart, &ct));
724: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, ct, 0, PETSC_DEFAULT, &feErr));
725: PetscCall(PetscObjectSetName((PetscObject)feErr, "Error"));
726: PetscCall(DMSetField(dmErr, 0, NULL, (PetscObject)feErr));
727: PetscCall(PetscFEDestroy(&feErr));
728: PetscCall(DMCreateDS(dmErr));
729: PetscCall(DMViewFromOptions(dmErr, NULL, "-dmerr_view"));
731: // Compute the flux norm
732: PetscCall(DMGetGlobalVector(dmErr, &fluxError));
733: PetscCall(PetscObjectSetName((PetscObject)fluxError, "Flux Error"));
734: PetscCall(DMGetGlobalVector(dmErr, &exactError));
735: PetscCall(PetscObjectSetName((PetscObject)exactError, "Analytical Error"));
736: PetscCall(DMPlexComputeL2FluxDiffVec(u, 0, mu, 0, fluxError));
737: {
738: PetscDS ds;
739: PetscSimplePointFn *func[2] = {NULL, NULL};
740: void *ctx[2] = {NULL, NULL};
742: PetscCall(DMGetDS(mdm, &ds));
743: PetscCall(PetscDSGetExactSolution(ds, 0, &func[0], &ctx[0]));
744: PetscCall(DMPlexComputeL2DiffVec(mdm, 0.0, func, ctx, mu, exactError));
745: }
746: PetscCall(VecNorm(fluxError, NORM_2, &fluxNorm));
747: PetscCall(VecNorm(exactError, NORM_2, &exactNorm));
748: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Flux error norm = %g\t Exact flux error norm = %g\n", (double)fluxNorm, (double)exactNorm));
749: PetscCall(VecViewFromOptions(fluxError, NULL, "-flux_error_vec_view"));
750: PetscCall(VecViewFromOptions(exactError, NULL, "-exact_error_vec_view"));
752: // Adaptive refinement based on calculated error
753: DM rdm;
754: VecTagger refineTag;
755: DMLabel adaptLabel;
756: IS refineIS;
757: Vec ref;
759: PetscCall(DMLabelCreate(PETSC_COMM_WORLD, "adapt", &adaptLabel));
760: PetscCall(DMLabelSetDefaultValue(adaptLabel, DM_ADAPT_COARSEN));
761: PetscCall(VecTaggerCreate(PETSC_COMM_WORLD, &refineTag));
762: PetscCall(VecTaggerSetFromOptions(refineTag));
763: PetscCall(VecTaggerSetUp(refineTag));
764: PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));
766: PetscCall(VecTaggerComputeIS(refineTag, fluxError, &refineIS, NULL));
767: PetscCall(VecTaggerDestroy(&refineTag));
768: PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
769: PetscCall(ISViewFromOptions(refineIS, NULL, "-refine_is_view"));
770: PetscCall(ISDestroy(&refineIS));
772: PetscCall(DMPlexCreateLabelField(dm, adaptLabel, &ref));
773: PetscCall(VecViewFromOptions(ref, NULL, "-refine_vec_view"));
774: PetscCall(VecDestroy(&ref));
776: // Mark adaptation phase with prefix ref_
777: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, "adapt_"));
778: PetscCall(DMAdaptLabel(dm, adaptLabel, &rdm));
779: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, NULL));
780: PetscCall(PetscObjectSetName((PetscObject)rdm, "Adaptively Refined DM"));
781: PetscCall(DMViewFromOptions(rdm, NULL, "-adapt_dm_view"));
782: PetscCall(DMDestroy(&rdm));
783: PetscCall(DMLabelDestroy(&adaptLabel));
785: // Destroy the error structures
786: PetscCall(DMRestoreGlobalVector(dmErr, &fluxError));
787: PetscCall(DMRestoreGlobalVector(dmErr, &exactError));
788: PetscCall(DMDestroy(&dmErr));
790: // Destroy the mixed structures
791: PetscCall(VecDestroy(&mu));
792: PetscCall(DMDestroy(&mdm));
793: PetscCall(SNESDestroy(&msnes));
794: }
796: // Destroy the primal structures
797: PetscCall(VecDestroy(&u));
798: PetscCall(DMDestroy(&dm));
799: PetscCall(SNESDestroy(&snes));
800: PetscCall(PetscBagDestroy(&user.param));
801: PetscCall(PetscFinalize());
802: return 0;
803: }
805: /*TEST
807: # Tests using the explicit code above
808: testset:
809: suffix: 2d_p2_rt0p0_byhand
810: requires: triangle
811: args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \
812: -primal_potential_petscspace_degree 2 \
813: -mixed_potential_petscdualspace_lagrange_continuity 0 \
814: -mixed_flux_petscspace_type ptrimmed \
815: -mixed_flux_petscspace_components 2 \
816: -mixed_flux_petscspace_ptrimmed_form_degree -1 \
817: -mixed_flux_petscdualspace_order 1 \
818: -mixed_flux_petscdualspace_form_degree -1 \
819: -mixed_flux_petscdualspace_lagrange_trimmed true \
820: -mixed_flux_petscfe_default_quadrature_order 2 \
821: -vec_tagger_type cdf -vec_tagger_box 0.9,1.0 \
822: -tag_view \
823: -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \
824: -dmsnes_check 0.001 -mixed_dmsnes_check 0.001 -pc_type jacobi -mixed_pc_type jacobi
825: test:
826: suffix: quadratic
827: args: -sol_type quadratic
828: test:
829: suffix: trig
830: args: -sol_type trig
831: test:
832: suffix: sensor
833: args: -sol_type sensor
835: # Tests using the embedded adaptor in SNES
836: testset:
837: suffix: 2d_p2_rt0p0
838: requires: triangle defined(PETSC_HAVE_EXECUTABLE_EXPORT)
839: args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \
840: -primal_potential_petscspace_degree 2 \
841: -mixed_potential_petscdualspace_lagrange_continuity 0 \
842: -mixed_flux_petscspace_type ptrimmed \
843: -mixed_flux_petscspace_components 2 \
844: -mixed_flux_petscspace_ptrimmed_form_degree -1 \
845: -mixed_flux_petscdualspace_order 1 \
846: -mixed_flux_petscdualspace_form_degree -1 \
847: -mixed_flux_petscdualspace_lagrange_trimmed true \
848: -mixed_flux_petscfe_default_quadrature_order 2 \
849: -by_hand 0 \
850: -refine_vec_tagger_type cdf -refine_vec_tagger_box 0.9,1.0 \
851: -snes_adapt_view \
852: -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \
853: -adaptor_criterion label -adaptor_type flux -adaptor_mixed_setup_function SetupMixed \
854: -snes_adapt_sequence 1 -pc_type jacobi -mixed_pc_type jacobi
855: test:
856: suffix: quadratic
857: args: -sol_type quadratic -adaptor_monitor_error
858: test:
859: suffix: trig
860: args: -sol_type trig -adaptor_monitor_error
861: test:
862: suffix: sensor
863: args: -sol_type sensor
865: # Tests using multiple adaptor loops
866: testset:
867: suffix: 2d_p2_rt0p0_a2
868: requires: triangle defined(PETSC_HAVE_EXECUTABLE_EXPORT)
869: args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \
870: -primal_potential_petscspace_degree 2 \
871: -mixed_potential_petscdualspace_lagrange_continuity 0 \
872: -mixed_flux_petscspace_type ptrimmed \
873: -mixed_flux_petscspace_components 2 \
874: -mixed_flux_petscspace_ptrimmed_form_degree -1 \
875: -mixed_flux_petscdualspace_order 1 \
876: -mixed_flux_petscdualspace_form_degree -1 \
877: -mixed_flux_petscdualspace_lagrange_trimmed true \
878: -mixed_flux_petscfe_default_quadrature_order 2 \
879: -by_hand 0 \
880: -num_adapt 2 \
881: -refine_vec_tagger_type cdf -refine_vec_tagger_box 0.9,1.0 \
882: -snes_adapt_view \
883: -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \
884: -adaptor_criterion label -adaptor_type gradient -adaptor_mixed_setup_function SetupMixed \
885: -snes_adapt_sequence 2 -pc_type jacobi \
886: -a1_refine_vec_tagger_type cdf -a1_refine_vec_tagger_box 0.9,1.0 \
887: -a1_snes_adapt_view \
888: -a1_adaptor_criterion label -a1_adaptor_type flux -a1_adaptor_mixed_setup_function SetupMixed \
889: -a1_snes_adapt_sequence 1 -a1_pc_type jacobi -a1_mixed_pc_type jacobi
890: test:
891: suffix: sensor
892: args: -sol_type sensor
894: TEST*/