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 <petscsnes.h>
91: #include <petscdmplex.h>
92: #include <petscdmadaptor.h>
93: #include <petscds.h>
94: #include <petscviewerhdf5.h>
95: #include <petscbag.h>
97: PETSC_EXTERN PetscErrorCode SetupMixed(DMAdaptor, DM);
99: typedef enum {
100: SOL_QUADRATIC,
101: SOL_TRIG,
102: SOL_SENSOR,
103: SOL_UNKNOWN,
104: NUM_SOL_TYPE
105: } SolType;
106: const char *SolTypeNames[NUM_SOL_TYPE + 4] = {"quadratic", "trig", "sensor", "unknown", "SolType", "SOL_", NULL};
108: typedef struct {
109: PetscBag param;
110: SolType solType;
111: PetscBool byHand;
112: PetscInt numAdapt;
113: } AppCtx;
115: typedef struct {
116: PetscReal k;
117: } Parameter;
119: /* Exact solution: u = x^2 + y^2 */
120: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
121: {
122: u[0] = 0.0;
123: for (PetscInt d = 0; d < dim; ++d) u[0] += x[d] * x[d];
124: return PETSC_SUCCESS;
125: }
126: /* Exact solution: q = (2x, 2y) */
127: static PetscErrorCode quadratic_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
128: {
129: for (PetscInt c = 0; c < Nc; ++c) u[c] = 2.0 * x[c];
130: return PETSC_SUCCESS;
131: }
133: /* Exact solution: u = sin( n \pi x ) * sin( n \pi y ) */
134: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
135: {
136: const PetscReal n = 5.5;
138: u[0] = 1.0;
139: for (PetscInt d = 0; d < dim; ++d) u[0] *= PetscSinReal(n * PETSC_PI * x[d]);
140: return PETSC_SUCCESS;
141: }
142: static PetscErrorCode trig_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
143: {
144: const PetscReal n = 5.5;
146: 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]);
147: return PETSC_SUCCESS;
148: }
150: /*
151: Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation:
153: f:[-1, 1]x[-1, 1] \to R,
154: f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy)
156: (mapped to have domain [0,1] x [0,1] in this case).
157: */
158: static PetscErrorCode sensor_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], PetscCtx ctx)
159: {
160: const PetscReal xref = 2. * x[0] - 1.;
161: const PetscReal yref = 2. * x[1] - 1.;
162: const PetscReal xy = xref * yref;
164: u[0] = PetscSinReal(50. * xy);
165: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) u[0] *= 0.01;
167: return PETSC_SUCCESS;
168: }
170: /* Flux is (cos(50xy) * 50y/100, cos(50xy) * 50x/100) if |xy| > 2\pi/50 else (cos(50xy) * 50y, cos(50xy) * 50x) */
171: static PetscErrorCode sensor_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], PetscCtx ctx)
172: {
173: const PetscReal xref = 2. * x[0] - 1.;
174: const PetscReal yref = 2. * x[1] - 1.;
175: const PetscReal xy = xref * yref;
177: u[0] = 50. * yref * PetscCosReal(50. * xy) * 2.0;
178: u[1] = 50. * xref * PetscCosReal(50. * xy) * 2.0;
179: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) {
180: u[0] *= 0.01;
181: u[1] *= 0.01;
182: }
183: return PETSC_SUCCESS;
184: }
186: /* We set up residuals and Jacobians for the primal problem. */
187: 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[])
188: {
189: f0[0] = 4.0;
190: }
192: 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[])
193: {
194: const PetscReal n = 5.5;
196: f0[0] = -2.0 * PetscSqr(n * PETSC_PI) * PetscSinReal(n * PETSC_PI * x[0]) * PetscSinReal(n * PETSC_PI * x[1]);
197: }
199: 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[])
200: {
201: const PetscReal xref = 2. * x[0] - 1.;
202: const PetscReal yref = 2. * x[1] - 1.;
203: const PetscReal xy = xref * yref;
205: f0[0] = -2500.0 * PetscSinReal(50. * xy) * (xref * xref + yref * yref) * 4.0;
206: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) f0[0] *= 0.01;
207: }
209: 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[])
210: {
211: const PetscReal k = PetscRealPart(constants[0]);
213: for (PetscInt d = 0; d < dim; ++d) f1[d] = k * u_x[uOff_x[0] + d];
214: }
216: 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[])
217: {
218: const PetscReal k = PetscRealPart(constants[0]);
219: PetscScalar flux;
221: PetscCallAbort(PETSC_COMM_SELF, quadratic_q(dim, t, x, dim, &flux, NULL));
222: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux * n[d];
223: }
225: 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[])
226: {
227: const PetscReal k = PetscRealPart(constants[0]);
228: PetscScalar flux;
230: PetscCallAbort(PETSC_COMM_SELF, trig_q(dim, t, x, dim, &flux, NULL));
231: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux * n[d];
232: }
234: 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[])
235: {
236: const PetscReal k = PetscRealPart(constants[0]);
237: PetscScalar flux[2];
239: PetscCallAbort(PETSC_COMM_SELF, sensor_q(dim, t, x, dim, flux, NULL));
240: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux[d] * n[d];
241: }
243: 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[])
244: {
245: const PetscReal k = PetscRealPart(constants[0]);
247: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = k;
248: }
250: /* Now we set up the residuals and Jacobians mixed problem. */
251: 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[])
252: {
253: f0[0] = 4.0;
254: for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d];
255: }
256: 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[])
257: {
258: PetscReal n = 5.5;
260: f0[0] = -2.0 * PetscSqr(n * PETSC_PI) * PetscSinReal(n * PETSC_PI * x[0]) * PetscSinReal(n * PETSC_PI * x[1]);
261: for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d];
262: }
263: 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[])
264: {
265: const PetscReal xref = 2. * x[0] - 1.;
266: const PetscReal yref = 2. * x[1] - 1.;
267: const PetscReal xy = xref * yref;
269: f0[0] = -2500.0 * PetscSinReal(50. * xy) * (xref * xref + yref * yref) * 4.0;
270: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) f0[0] *= 0.01;
271: for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d];
272: }
274: 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[])
275: {
276: for (PetscInt d = 0; d < dim; d++) f0[d] = u[uOff[0] + d];
277: }
279: 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[])
280: {
281: const PetscReal k = PetscRealPart(constants[0]);
283: for (PetscInt d = 0; d < dim; d++) f1[d * dim + d] = k * u[uOff[1]];
284: }
286: 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[])
287: {
288: const PetscReal k = PetscRealPart(constants[0]);
289: PetscScalar potential;
291: PetscCallAbort(PETSC_COMM_SELF, quadratic_u(dim, t, x, dim, &potential, NULL));
292: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * potential * n[d];
293: }
295: 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[])
296: {
297: const PetscReal k = PetscRealPart(constants[0]);
298: PetscScalar potential;
300: PetscCallAbort(PETSC_COMM_SELF, trig_u(dim, t, x, dim, &potential, NULL));
301: for (PetscInt d = 0; d < dim; ++d) f0[d * dim + d] = -k * potential * n[d];
302: }
304: 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[])
305: {
306: const PetscReal k = PetscRealPart(constants[0]);
307: PetscScalar potential;
309: PetscCallAbort(PETSC_COMM_SELF, sensor_u(dim, t, x, dim, &potential, NULL));
310: for (PetscInt d = 0; d < dim; ++d) f0[d * dim + d] = -k * potential * n[d];
311: }
313: /* <v, \nabla\cdot q> */
314: 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[])
315: {
316: for (PetscInt d = 0; d < dim; d++) g1[d * dim + d] = -1.0;
317: }
319: /* < t, q> */
320: 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[])
321: {
322: for (PetscInt d = 0; d < dim; d++) g0[d * dim + d] = 1.0;
323: }
325: /* <\nabla\cdot t, u> = <\nabla t, Iu> */
326: 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[])
327: {
328: const PetscReal k = PetscRealPart(constants[0]);
330: for (PetscInt d = 0; d < dim; d++) g2[d * dim + d] = k;
331: }
333: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *user)
334: {
335: PetscFunctionBeginUser;
336: PetscOptionsBegin(comm, "", "Flux norm error in primal poisson problem Options", "DMPLEX");
337: user->byHand = PETSC_TRUE;
338: user->numAdapt = 1;
339: user->solType = SOL_QUADRATIC;
341: PetscCall(PetscOptionsGetBool(NULL, NULL, "-by_hand", &user->byHand, NULL));
342: PetscCall(PetscOptionsGetInt(NULL, NULL, "-num_adapt", &user->numAdapt, NULL));
343: PetscCall(PetscOptionsEnum("-sol_type", "Type of exact solution", "ex27.c", SolTypeNames, (PetscEnum)user->solType, (PetscEnum *)&user->solType, NULL));
344: PetscOptionsEnd();
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: static PetscErrorCode SetupParameters(PetscBag bag, AppCtx *user)
349: {
350: Parameter *param;
352: PetscFunctionBeginUser;
353: PetscCall(PetscBagGetData(bag, ¶m));
354: PetscCall(PetscBagSetName(bag, "par", "Poisson parameters"));
355: PetscCall(PetscBagRegisterReal(bag, ¶m->k, 1.0, "k", "Thermal conductivity"));
356: PetscCall(PetscBagSetFromOptions(bag));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
361: {
362: PetscFunctionBeginUser;
363: PetscCall(DMCreate(comm, dm));
364: PetscCall(DMSetType(*dm, DMPLEX));
365: PetscCall(DMSetFromOptions(*dm));
366: PetscCall(DMSetApplicationContext(*dm, &user));
367: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
372: {
373: PetscDS ds;
374: DMLabel label;
375: PetscInt id, bd;
376: Parameter *param;
377: PetscWeakForm wf;
379: PetscFunctionBeginUser;
380: PetscCall(DMGetDS(dm, &ds));
381: PetscCall(DMGetLabel(dm, "marker", &label));
383: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_primal_uu));
385: switch (user->solType) {
386: case SOL_QUADRATIC:
387: PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_primal, f1_primal));
389: id = 1;
390: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_u, NULL, user, NULL));
392: id = 2;
393: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
394: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
395: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_primal, 0, NULL));
397: id = 3;
398: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_u, NULL, user, NULL));
400: id = 4;
401: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
402: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
403: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_primal, 0, NULL));
405: PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_u, user));
406: break;
407: case SOL_TRIG:
408: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_primal, f1_primal));
410: id = 1;
411: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL));
413: id = 2;
414: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
415: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
416: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_primal, 0, NULL));
418: id = 3;
419: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL));
421: id = 4;
422: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
423: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
424: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_primal, 0, NULL));
426: PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user));
427: break;
428: case SOL_SENSOR:
429: PetscCall(PetscDSSetResidual(ds, 0, f0_sensor_primal, f1_primal));
431: id = 1;
432: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)sensor_u, NULL, user, NULL));
434: id = 2;
435: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
436: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
437: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_primal, 0, NULL));
439: id = 3;
440: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)sensor_u, NULL, user, NULL));
442: id = 4;
443: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
444: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
445: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_primal, 0, NULL));
447: PetscCall(PetscDSSetExactSolution(ds, 0, sensor_u, user));
448: break;
449: default:
450: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]);
451: }
453: /* Setup constants */
454: {
455: PetscCall(PetscBagGetData(user->param, ¶m));
456: PetscScalar constants[1];
458: constants[0] = param->k;
460: PetscCall(PetscDSSetConstants(ds, 1, constants));
461: }
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: static PetscErrorCode SetupPrimalDiscretization(DM dm, AppCtx *user)
466: {
467: DM cdm = dm;
468: PetscFE fe[1];
469: DMPolytopeType ct;
470: PetscInt dim, cStart;
471: MPI_Comm comm;
473: PetscFunctionBeginUser;
474: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
475: PetscCall(DMGetDimension(dm, &dim));
477: /* Create finite element */
478: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
479: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
480: PetscCall(PetscFECreateByCell(comm, dim, 1, ct, "primal_potential_", PETSC_DEFAULT, &fe[0]));
481: PetscCall(PetscObjectSetName((PetscObject)fe[0], "primal potential"));
483: /* Set discretization and boundary conditions for each mesh */
484: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
485: PetscCall(DMCreateDS(dm));
486: while (cdm) {
487: PetscCall(DMCopyDisc(dm, cdm));
488: PetscCall(DMGetCoarseDM(cdm, &cdm));
489: }
491: PetscCall(PetscFEDestroy(&fe[0]));
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: static PetscErrorCode SetupMixedProblem(DM dm, AppCtx *user)
496: {
497: PetscDS ds;
498: DMLabel label;
499: PetscInt id, bd;
500: Parameter *param;
501: PetscWeakForm wf;
503: PetscFunctionBeginUser;
504: PetscCall(DMGetDS(dm, &ds));
505: PetscCall(DMGetLabel(dm, "marker", &label));
507: /* Residual terms */
508: PetscCall(PetscDSSetResidual(ds, 0, f0_mixed_q, f1_mixed_q));
510: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_mixed_qq, NULL, NULL, NULL));
511: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_mixed_qu, NULL));
513: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_mixed_uq, NULL, NULL));
515: switch (user->solType) {
516: case SOL_QUADRATIC:
517: PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_quadratic_u, NULL));
519: id = 1;
520: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
521: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
522: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_mixed_q, 0, NULL));
524: id = 2;
525: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_q, NULL, user, NULL));
527: id = 4;
528: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_q, NULL, user, NULL));
530: id = 3;
531: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
532: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
533: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_mixed_q, 0, NULL));
535: PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_q, user));
536: PetscCall(PetscDSSetExactSolution(ds, 1, quadratic_u, user));
537: break;
538: case SOL_TRIG:
539: PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_trig_u, NULL));
541: id = 1;
542: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
543: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
544: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_mixed_q, 0, NULL));
546: id = 2;
547: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)trig_q, NULL, user, NULL));
549: id = 4;
550: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)trig_q, NULL, user, NULL));
552: id = 3;
553: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
554: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
555: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_mixed_q, 0, NULL));
557: PetscCall(PetscDSSetExactSolution(ds, 0, trig_q, user));
558: PetscCall(PetscDSSetExactSolution(ds, 1, trig_u, user));
559: break;
560: case SOL_SENSOR:
561: PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_sensor_u, NULL));
563: id = 1;
564: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
565: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
566: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_mixed_q, 0, NULL));
568: id = 2;
569: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)sensor_q, NULL, user, NULL));
571: id = 4;
572: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)sensor_q, NULL, user, NULL));
574: id = 3;
575: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
576: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
577: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_mixed_q, 0, NULL));
579: PetscCall(PetscDSSetExactSolution(ds, 0, sensor_q, user));
580: PetscCall(PetscDSSetExactSolution(ds, 1, sensor_u, user));
581: break;
582: default:
583: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]);
584: }
586: /* Setup constants */
587: {
588: PetscCall(PetscBagGetData(user->param, ¶m));
589: PetscScalar constants[1];
591: constants[0] = param->k;
593: PetscCall(PetscDSSetConstants(ds, 1, constants));
594: }
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: static PetscErrorCode SetupMixedDiscretization(DM dm, AppCtx *user)
599: {
600: DM cdm = dm;
601: PetscFE fe[2];
602: DMPolytopeType ct;
603: PetscInt dim, cStart;
604: MPI_Comm comm;
606: PetscFunctionBeginUser;
607: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
608: PetscCall(DMGetDimension(dm, &dim));
610: /* Create finite element */
611: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
612: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
613: PetscCall(PetscFECreateByCell(comm, dim, dim, ct, "mixed_flux_", PETSC_DEFAULT, &fe[0]));
614: PetscCall(PetscObjectSetName((PetscObject)fe[0], "mixed flux"));
615: /* NOTE: Set the same quadrature order as the primal problem here or use the command line option. */
617: PetscCall(PetscFECreateByCell(comm, dim, 1, ct, "mixed_potential_", PETSC_DEFAULT, &fe[1]));
618: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
619: PetscCall(PetscObjectSetName((PetscObject)fe[1], "mixed potential"));
621: /* Set discretization and boundary conditions for each mesh */
622: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
623: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
624: PetscCall(DMCreateDS(dm));
625: while (cdm) {
626: PetscCall(DMCopyDisc(dm, cdm));
627: PetscCall(DMGetCoarseDM(cdm, &cdm));
628: }
629: PetscCall(PetscFEDestroy(&fe[0]));
630: PetscCall(PetscFEDestroy(&fe[1]));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: PetscErrorCode SetupMixed(DMAdaptor adaptor, DM mdm)
635: {
636: AppCtx *ctx;
638: PetscFunctionBeginUser;
639: PetscCall(DMGetApplicationContext(mdm, &ctx));
640: PetscCall(SetupMixedDiscretization(mdm, ctx));
641: PetscCall(SetupMixedProblem(mdm, ctx));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: int main(int argc, char **argv)
646: {
647: DM dm, mdm; /* problem specification */
648: SNES snes, msnes; /* nonlinear solvers */
649: Vec u, mu; /* solution vectors */
650: Vec fluxError, exactError; /* Element wise error vector */
651: PetscReal fluxNorm, exactNorm; /* Flux error norm */
652: AppCtx user; /* user-defined work context */
654: PetscFunctionBeginUser;
655: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
656: PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.param));
657: PetscCall(SetupParameters(user.param, &user));
658: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
660: // Set up and solve primal problem
661: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
662: PetscCall(DMSetApplicationContext(dm, &user));
663: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
664: PetscCall(SNESSetDM(snes, dm));
666: PetscCall(SetupPrimalDiscretization(dm, &user));
667: PetscCall(SetupPrimalProblem(dm, &user));
668: PetscCall(DMCreateGlobalVector(dm, &u));
669: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
670: PetscCall(SNESSetFromOptions(snes));
671: PetscCall(DMSNESCheckFromOptions(snes, u));
673: for (PetscInt a = 0; a < user.numAdapt; ++a) {
674: if (a > 0) {
675: char prefix[16];
677: PetscCall(PetscSNPrintf(prefix, 16, "a%d_", (int)a));
678: PetscCall(SNESSetOptionsPrefix(snes, prefix));
679: }
680: PetscCall(SNESSolve(snes, NULL, u));
682: // Needed if you allow SNES to refine
683: PetscCall(SNESGetSolution(snes, &u));
684: PetscCall(VecGetDM(u, &dm));
685: }
687: PetscCall(PetscObjectSetName((PetscObject)u, "Primal Solution "));
688: PetscCall(VecViewFromOptions(u, NULL, "-primal_sol_vec_view"));
690: if (user.byHand) {
691: // Set up and solve mixed problem
692: PetscCall(DMClone(dm, &mdm));
693: PetscCall(SNESCreate(PETSC_COMM_WORLD, &msnes));
694: PetscCall(SNESSetOptionsPrefix(msnes, "mixed_"));
695: PetscCall(SNESSetDM(msnes, mdm));
697: PetscCall(SetupMixedDiscretization(mdm, &user));
698: PetscCall(SetupMixedProblem(mdm, &user));
699: PetscCall(DMCreateGlobalVector(mdm, &mu));
700: PetscCall(DMPlexSetSNESLocalFEM(mdm, PETSC_FALSE, &user));
701: PetscCall(SNESSetFromOptions(msnes));
703: PetscCall(DMSNESCheckFromOptions(msnes, mu));
704: PetscCall(SNESSolve(msnes, NULL, mu));
705: PetscCall(PetscObjectSetName((PetscObject)mu, "Mixed Solution "));
706: PetscCall(VecViewFromOptions(mu, NULL, "-mixed_sol_vec_view"));
708: // Create the error space of piecewise constants
709: DM dmErr;
710: PetscFE feErr;
711: DMPolytopeType ct;
712: PetscInt dim, cStart;
714: PetscCall(DMClone(dm, &dmErr));
715: PetscCall(DMGetDimension(dmErr, &dim));
716: PetscCall(DMPlexGetHeightStratum(dmErr, 0, &cStart, NULL));
717: PetscCall(DMPlexGetCellType(dmErr, cStart, &ct));
718: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, ct, 0, PETSC_DEFAULT, &feErr));
719: PetscCall(PetscObjectSetName((PetscObject)feErr, "Error"));
720: PetscCall(DMSetField(dmErr, 0, NULL, (PetscObject)feErr));
721: PetscCall(PetscFEDestroy(&feErr));
722: PetscCall(DMCreateDS(dmErr));
723: PetscCall(DMViewFromOptions(dmErr, NULL, "-dmerr_view"));
725: // Compute the flux norm
726: PetscCall(DMGetGlobalVector(dmErr, &fluxError));
727: PetscCall(PetscObjectSetName((PetscObject)fluxError, "Flux Error"));
728: PetscCall(DMGetGlobalVector(dmErr, &exactError));
729: PetscCall(PetscObjectSetName((PetscObject)exactError, "Analytical Error"));
730: PetscCall(DMPlexComputeL2FluxDiffVec(u, 0, mu, 0, fluxError));
731: {
732: PetscDS ds;
733: PetscSimplePointFn *func[2] = {NULL, NULL};
734: void *ctx[2] = {NULL, NULL};
736: PetscCall(DMGetDS(mdm, &ds));
737: PetscCall(PetscDSGetExactSolution(ds, 0, &func[0], &ctx[0]));
738: PetscCall(DMPlexComputeL2DiffVec(mdm, 0.0, func, ctx, mu, exactError));
739: }
740: PetscCall(VecNorm(fluxError, NORM_2, &fluxNorm));
741: PetscCall(VecNorm(exactError, NORM_2, &exactNorm));
742: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Flux error norm = %g\t Exact flux error norm = %g\n", (double)fluxNorm, (double)exactNorm));
743: PetscCall(VecViewFromOptions(fluxError, NULL, "-flux_error_vec_view"));
744: PetscCall(VecViewFromOptions(exactError, NULL, "-exact_error_vec_view"));
746: // Adaptive refinement based on calculated error
747: DM rdm;
748: VecTagger refineTag;
749: DMLabel adaptLabel;
750: IS refineIS;
751: Vec ref;
753: PetscCall(DMLabelCreate(PETSC_COMM_WORLD, "adapt", &adaptLabel));
754: PetscCall(DMLabelSetDefaultValue(adaptLabel, DM_ADAPT_COARSEN));
755: PetscCall(VecTaggerCreate(PETSC_COMM_WORLD, &refineTag));
756: PetscCall(VecTaggerSetFromOptions(refineTag));
757: PetscCall(VecTaggerSetUp(refineTag));
758: PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));
760: PetscCall(VecTaggerComputeIS(refineTag, fluxError, &refineIS, NULL));
761: PetscCall(VecTaggerDestroy(&refineTag));
762: PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
763: PetscCall(ISViewFromOptions(refineIS, NULL, "-refine_is_view"));
764: PetscCall(ISDestroy(&refineIS));
766: PetscCall(DMPlexCreateLabelField(dm, adaptLabel, &ref));
767: PetscCall(VecViewFromOptions(ref, NULL, "-refine_vec_view"));
768: PetscCall(VecDestroy(&ref));
770: // Mark adaptation phase with prefix ref_
771: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, "adapt_"));
772: PetscCall(DMAdaptLabel(dm, adaptLabel, &rdm));
773: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, NULL));
774: PetscCall(PetscObjectSetName((PetscObject)rdm, "Adaptively Refined DM"));
775: PetscCall(DMViewFromOptions(rdm, NULL, "-adapt_dm_view"));
776: PetscCall(DMDestroy(&rdm));
777: PetscCall(DMLabelDestroy(&adaptLabel));
779: // Destroy the error structures
780: PetscCall(DMRestoreGlobalVector(dmErr, &fluxError));
781: PetscCall(DMRestoreGlobalVector(dmErr, &exactError));
782: PetscCall(DMDestroy(&dmErr));
784: // Destroy the mixed structures
785: PetscCall(VecDestroy(&mu));
786: PetscCall(DMDestroy(&mdm));
787: PetscCall(SNESDestroy(&msnes));
788: }
790: // Destroy the primal structures
791: PetscCall(VecDestroy(&u));
792: PetscCall(DMDestroy(&dm));
793: PetscCall(SNESDestroy(&snes));
794: PetscCall(PetscBagDestroy(&user.param));
795: PetscCall(PetscFinalize());
796: return 0;
797: }
799: /*TEST
801: # Tests using the explicit code above
802: testset:
803: suffix: 2d_p2_rt0p0_byhand
804: requires: triangle
805: args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \
806: -primal_potential_petscspace_degree 2 \
807: -mixed_potential_petscdualspace_lagrange_continuity 0 \
808: -mixed_flux_petscspace_type ptrimmed \
809: -mixed_flux_petscspace_components 2 \
810: -mixed_flux_petscspace_ptrimmed_form_degree -1 \
811: -mixed_flux_petscdualspace_order 1 \
812: -mixed_flux_petscdualspace_form_degree -1 \
813: -mixed_flux_petscdualspace_lagrange_trimmed true \
814: -mixed_flux_petscfe_default_quadrature_order 2 \
815: -vec_tagger_type cdf -vec_tagger_box 0.9,1.0 \
816: -tag_view \
817: -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \
818: -dmsnes_check 0.001 -mixed_dmsnes_check 0.001 -pc_type jacobi -mixed_pc_type jacobi
819: test:
820: suffix: quadratic
821: args: -sol_type quadratic
822: test:
823: suffix: trig
824: args: -sol_type trig
825: test:
826: suffix: sensor
827: args: -sol_type sensor
829: # Tests using the embedded adaptor in SNES
830: testset:
831: suffix: 2d_p2_rt0p0
832: requires: triangle defined(PETSC_HAVE_EXECUTABLE_EXPORT)
833: args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \
834: -primal_potential_petscspace_degree 2 \
835: -mixed_potential_petscdualspace_lagrange_continuity 0 \
836: -mixed_flux_petscspace_type ptrimmed \
837: -mixed_flux_petscspace_components 2 \
838: -mixed_flux_petscspace_ptrimmed_form_degree -1 \
839: -mixed_flux_petscdualspace_order 1 \
840: -mixed_flux_petscdualspace_form_degree -1 \
841: -mixed_flux_petscdualspace_lagrange_trimmed true \
842: -mixed_flux_petscfe_default_quadrature_order 2 \
843: -by_hand 0 \
844: -refine_vec_tagger_type cdf -refine_vec_tagger_box 0.9,1.0 \
845: -snes_adapt_view \
846: -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \
847: -adaptor_criterion label -adaptor_type flux -adaptor_mixed_setup_function SetupMixed \
848: -snes_adapt_sequence 1 -pc_type jacobi -mixed_pc_type jacobi
849: test:
850: suffix: quadratic
851: args: -sol_type quadratic -adaptor_monitor_error
852: test:
853: suffix: trig
854: args: -sol_type trig -adaptor_monitor_error
855: test:
856: suffix: sensor
857: args: -sol_type sensor
859: # Tests using multiple adaptor loops
860: testset:
861: suffix: 2d_p2_rt0p0_a2
862: requires: triangle defined(PETSC_HAVE_EXECUTABLE_EXPORT)
863: args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \
864: -primal_potential_petscspace_degree 2 \
865: -mixed_potential_petscdualspace_lagrange_continuity 0 \
866: -mixed_flux_petscspace_type ptrimmed \
867: -mixed_flux_petscspace_components 2 \
868: -mixed_flux_petscspace_ptrimmed_form_degree -1 \
869: -mixed_flux_petscdualspace_order 1 \
870: -mixed_flux_petscdualspace_form_degree -1 \
871: -mixed_flux_petscdualspace_lagrange_trimmed true \
872: -mixed_flux_petscfe_default_quadrature_order 2 \
873: -by_hand 0 \
874: -num_adapt 2 \
875: -refine_vec_tagger_type cdf -refine_vec_tagger_box 0.9,1.0 \
876: -snes_adapt_view \
877: -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \
878: -adaptor_criterion label -adaptor_type gradient -adaptor_mixed_setup_function SetupMixed \
879: -snes_adapt_sequence 2 -pc_type jacobi \
880: -a1_refine_vec_tagger_type cdf -a1_refine_vec_tagger_box 0.9,1.0 \
881: -a1_snes_adapt_view \
882: -a1_adaptor_criterion label -a1_adaptor_type flux -a1_adaptor_mixed_setup_function SetupMixed \
883: -a1_snes_adapt_sequence 1 -a1_pc_type jacobi -a1_mixed_pc_type jacobi
884: test:
885: suffix: sensor
886: args: -sol_type sensor
888: TEST*/