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, PetscCtx ctx)
122: {
123: u[0] = 0.0;
124: for (PetscInt d = 0; d < dim; ++d) u[0] += x[d] * x[d];
125: return PETSC_SUCCESS;
126: }
127: /* Exact solution: q = (2x, 2y) */
128: static PetscErrorCode quadratic_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
129: {
130: for (PetscInt c = 0; c < Nc; ++c) u[c] = 2.0 * x[c];
131: return PETSC_SUCCESS;
132: }
134: /* Exact solution: u = sin( n \pi x ) * sin( n \pi y ) */
135: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
136: {
137: const PetscReal n = 5.5;
139: u[0] = 1.0;
140: for (PetscInt d = 0; d < dim; ++d) u[0] *= PetscSinReal(n * PETSC_PI * x[d]);
141: return PETSC_SUCCESS;
142: }
143: static PetscErrorCode trig_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
144: {
145: const PetscReal n = 5.5;
147: 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]);
148: return PETSC_SUCCESS;
149: }
151: /*
152: Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation:
154: f:[-1, 1]x[-1, 1] \to R,
155: f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy)
157: (mapped to have domain [0,1] x [0,1] in this case).
158: */
159: static PetscErrorCode sensor_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], PetscCtx ctx)
160: {
161: const PetscReal xref = 2. * x[0] - 1.;
162: const PetscReal yref = 2. * x[1] - 1.;
163: const PetscReal xy = xref * yref;
165: u[0] = PetscSinReal(50. * xy);
166: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) u[0] *= 0.01;
168: return PETSC_SUCCESS;
169: }
171: /* Flux is (cos(50xy) * 50y/100, cos(50xy) * 50x/100) if |xy| > 2\pi/50 else (cos(50xy) * 50y, cos(50xy) * 50x) */
172: static PetscErrorCode sensor_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], PetscCtx ctx)
173: {
174: const PetscReal xref = 2. * x[0] - 1.;
175: const PetscReal yref = 2. * x[1] - 1.;
176: const PetscReal xy = xref * yref;
178: u[0] = 50. * yref * PetscCosReal(50. * xy) * 2.0;
179: u[1] = 50. * xref * PetscCosReal(50. * xy) * 2.0;
180: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) {
181: u[0] *= 0.01;
182: u[1] *= 0.01;
183: }
184: return PETSC_SUCCESS;
185: }
187: /* We set up residuals and Jacobians for the primal problem. */
188: 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[])
189: {
190: f0[0] = 4.0;
191: }
193: 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[])
194: {
195: const PetscReal n = 5.5;
197: f0[0] = -2.0 * PetscSqr(n * PETSC_PI) * PetscSinReal(n * PETSC_PI * x[0]) * PetscSinReal(n * PETSC_PI * x[1]);
198: }
200: 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[])
201: {
202: const PetscReal xref = 2. * x[0] - 1.;
203: const PetscReal yref = 2. * x[1] - 1.;
204: const PetscReal xy = xref * yref;
206: f0[0] = -2500.0 * PetscSinReal(50. * xy) * (xref * xref + yref * yref) * 4.0;
207: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) f0[0] *= 0.01;
208: }
210: 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[])
211: {
212: const PetscReal k = PetscRealPart(constants[0]);
214: for (PetscInt d = 0; d < dim; ++d) f1[d] = k * u_x[uOff_x[0] + d];
215: }
217: 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[])
218: {
219: const PetscReal k = PetscRealPart(constants[0]);
220: PetscScalar flux;
222: PetscCallAbort(PETSC_COMM_SELF, quadratic_q(dim, t, x, dim, &flux, NULL));
223: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux * n[d];
224: }
226: 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[])
227: {
228: const PetscReal k = PetscRealPart(constants[0]);
229: PetscScalar flux;
231: PetscCallAbort(PETSC_COMM_SELF, trig_q(dim, t, x, dim, &flux, NULL));
232: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux * n[d];
233: }
235: 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[])
236: {
237: const PetscReal k = PetscRealPart(constants[0]);
238: PetscScalar flux[2];
240: PetscCallAbort(PETSC_COMM_SELF, sensor_q(dim, t, x, dim, flux, NULL));
241: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux[d] * n[d];
242: }
244: 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[])
245: {
246: const PetscReal k = PetscRealPart(constants[0]);
248: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = k;
249: }
251: /* Now we set up the residuals and Jacobians mixed problem. */
252: 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[])
253: {
254: f0[0] = 4.0;
255: for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d];
256: }
257: 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[])
258: {
259: PetscReal n = 5.5;
261: f0[0] = -2.0 * PetscSqr(n * PETSC_PI) * PetscSinReal(n * PETSC_PI * x[0]) * PetscSinReal(n * PETSC_PI * x[1]);
262: for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d];
263: }
264: 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[])
265: {
266: const PetscReal xref = 2. * x[0] - 1.;
267: const PetscReal yref = 2. * x[1] - 1.;
268: const PetscReal xy = xref * yref;
270: f0[0] = -2500.0 * PetscSinReal(50. * xy) * (xref * xref + yref * yref) * 4.0;
271: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) f0[0] *= 0.01;
272: for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d];
273: }
275: 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[])
276: {
277: for (PetscInt d = 0; d < dim; d++) f0[d] = u[uOff[0] + d];
278: }
280: 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[])
281: {
282: const PetscReal k = PetscRealPart(constants[0]);
284: for (PetscInt d = 0; d < dim; d++) f1[d * dim + d] = k * u[uOff[1]];
285: }
287: 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[])
288: {
289: const PetscReal k = PetscRealPart(constants[0]);
290: PetscScalar potential;
292: PetscCallAbort(PETSC_COMM_SELF, quadratic_u(dim, t, x, dim, &potential, NULL));
293: for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * potential * n[d];
294: }
296: 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[])
297: {
298: const PetscReal k = PetscRealPart(constants[0]);
299: PetscScalar potential;
301: PetscCallAbort(PETSC_COMM_SELF, trig_u(dim, t, x, dim, &potential, NULL));
302: for (PetscInt d = 0; d < dim; ++d) f0[d * dim + d] = -k * potential * n[d];
303: }
305: 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[])
306: {
307: const PetscReal k = PetscRealPart(constants[0]);
308: PetscScalar potential;
310: PetscCallAbort(PETSC_COMM_SELF, sensor_u(dim, t, x, dim, &potential, NULL));
311: for (PetscInt d = 0; d < dim; ++d) f0[d * dim + d] = -k * potential * n[d];
312: }
314: /* <v, \nabla\cdot q> */
315: 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[])
316: {
317: for (PetscInt d = 0; d < dim; d++) g1[d * dim + d] = -1.0;
318: }
320: /* < t, q> */
321: 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[])
322: {
323: for (PetscInt d = 0; d < dim; d++) g0[d * dim + d] = 1.0;
324: }
326: /* <\nabla\cdot t, u> = <\nabla t, Iu> */
327: 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[])
328: {
329: const PetscReal k = PetscRealPart(constants[0]);
331: for (PetscInt d = 0; d < dim; d++) g2[d * dim + d] = k;
332: }
334: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *user)
335: {
336: PetscFunctionBeginUser;
337: PetscOptionsBegin(comm, "", "Flux norm error in primal poisson problem Options", "DMPLEX");
338: user->byHand = PETSC_TRUE;
339: user->numAdapt = 1;
340: user->solType = SOL_QUADRATIC;
342: PetscCall(PetscOptionsGetBool(NULL, NULL, "-by_hand", &user->byHand, NULL));
343: PetscCall(PetscOptionsGetInt(NULL, NULL, "-num_adapt", &user->numAdapt, NULL));
344: PetscCall(PetscOptionsEnum("-sol_type", "Type of exact solution", "ex27.c", SolTypeNames, (PetscEnum)user->solType, (PetscEnum *)&user->solType, NULL));
345: PetscOptionsEnd();
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: static PetscErrorCode SetupParameters(PetscBag bag, AppCtx *user)
350: {
351: Parameter *param;
353: PetscFunctionBeginUser;
354: PetscCall(PetscBagGetData(bag, ¶m));
355: PetscCall(PetscBagSetName(bag, "par", "Poisson parameters"));
356: PetscCall(PetscBagRegisterReal(bag, ¶m->k, 1.0, "k", "Thermal conductivity"));
357: PetscCall(PetscBagSetFromOptions(bag));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
362: {
363: PetscFunctionBeginUser;
364: PetscCall(DMCreate(comm, dm));
365: PetscCall(DMSetType(*dm, DMPLEX));
366: PetscCall(DMSetFromOptions(*dm));
367: PetscCall(DMSetApplicationContext(*dm, &user));
368: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
373: {
374: PetscDS ds;
375: DMLabel label;
376: PetscInt id, bd;
377: Parameter *param;
378: PetscWeakForm wf;
380: PetscFunctionBeginUser;
381: PetscCall(DMGetDS(dm, &ds));
382: PetscCall(DMGetLabel(dm, "marker", &label));
384: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_primal_uu));
386: switch (user->solType) {
387: case SOL_QUADRATIC:
388: PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_primal, f1_primal));
390: id = 1;
391: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_u, NULL, user, NULL));
393: id = 2;
394: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
395: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
396: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_primal, 0, NULL));
398: id = 3;
399: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_u, NULL, user, NULL));
401: id = 4;
402: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
403: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
404: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_primal, 0, NULL));
406: PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_u, user));
407: break;
408: case SOL_TRIG:
409: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_primal, f1_primal));
411: id = 1;
412: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL));
414: id = 2;
415: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
416: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
417: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_primal, 0, NULL));
419: id = 3;
420: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL));
422: id = 4;
423: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
424: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
425: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_primal, 0, NULL));
427: PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user));
428: break;
429: case SOL_SENSOR:
430: PetscCall(PetscDSSetResidual(ds, 0, f0_sensor_primal, f1_primal));
432: id = 1;
433: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)sensor_u, NULL, user, NULL));
435: id = 2;
436: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
437: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
438: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_primal, 0, NULL));
440: id = 3;
441: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)sensor_u, NULL, user, NULL));
443: id = 4;
444: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
445: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
446: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_primal, 0, NULL));
448: PetscCall(PetscDSSetExactSolution(ds, 0, sensor_u, user));
449: break;
450: default:
451: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]);
452: }
454: /* Setup constants */
455: {
456: PetscCall(PetscBagGetData(user->param, ¶m));
457: PetscScalar constants[1];
459: constants[0] = param->k;
461: PetscCall(PetscDSSetConstants(ds, 1, constants));
462: }
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: static PetscErrorCode SetupPrimalDiscretization(DM dm, AppCtx *user)
467: {
468: DM cdm = dm;
469: PetscFE fe[1];
470: DMPolytopeType ct;
471: PetscInt dim, cStart;
472: MPI_Comm comm;
474: PetscFunctionBeginUser;
475: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
476: PetscCall(DMGetDimension(dm, &dim));
478: /* Create finite element */
479: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
480: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
481: PetscCall(PetscFECreateByCell(comm, dim, 1, ct, "primal_potential_", PETSC_DEFAULT, &fe[0]));
482: PetscCall(PetscObjectSetName((PetscObject)fe[0], "primal potential"));
484: /* Set discretization and boundary conditions for each mesh */
485: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
486: PetscCall(DMCreateDS(dm));
487: while (cdm) {
488: PetscCall(DMCopyDisc(dm, cdm));
489: PetscCall(DMGetCoarseDM(cdm, &cdm));
490: }
492: PetscCall(PetscFEDestroy(&fe[0]));
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: static PetscErrorCode SetupMixedProblem(DM dm, AppCtx *user)
497: {
498: PetscDS ds;
499: DMLabel label;
500: PetscInt id, bd;
501: Parameter *param;
502: PetscWeakForm wf;
504: PetscFunctionBeginUser;
505: PetscCall(DMGetDS(dm, &ds));
506: PetscCall(DMGetLabel(dm, "marker", &label));
508: /* Residual terms */
509: PetscCall(PetscDSSetResidual(ds, 0, f0_mixed_q, f1_mixed_q));
511: PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_mixed_qq, NULL, NULL, NULL));
512: PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_mixed_qu, NULL));
514: PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_mixed_uq, NULL, NULL));
516: switch (user->solType) {
517: case SOL_QUADRATIC:
518: PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_quadratic_u, NULL));
520: id = 1;
521: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
522: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
523: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_mixed_q, 0, NULL));
525: id = 2;
526: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_q, NULL, user, NULL));
528: id = 4;
529: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_q, NULL, user, NULL));
531: id = 3;
532: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
533: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
534: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_mixed_q, 0, NULL));
536: PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_q, user));
537: PetscCall(PetscDSSetExactSolution(ds, 1, quadratic_u, user));
538: break;
539: case SOL_TRIG:
540: PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_trig_u, NULL));
542: id = 1;
543: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
544: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
545: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_mixed_q, 0, NULL));
547: id = 2;
548: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)trig_q, NULL, user, NULL));
550: id = 4;
551: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)trig_q, NULL, user, NULL));
553: id = 3;
554: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
555: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
556: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_mixed_q, 0, NULL));
558: PetscCall(PetscDSSetExactSolution(ds, 0, trig_q, user));
559: PetscCall(PetscDSSetExactSolution(ds, 1, trig_u, user));
560: break;
561: case SOL_SENSOR:
562: PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_sensor_u, NULL));
564: id = 1;
565: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
566: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
567: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_mixed_q, 0, NULL));
569: id = 2;
570: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)sensor_q, NULL, user, NULL));
572: id = 4;
573: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)sensor_q, NULL, user, NULL));
575: id = 3;
576: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd));
577: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
578: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_mixed_q, 0, NULL));
580: PetscCall(PetscDSSetExactSolution(ds, 0, sensor_q, user));
581: PetscCall(PetscDSSetExactSolution(ds, 1, sensor_u, user));
582: break;
583: default:
584: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]);
585: }
587: /* Setup constants */
588: {
589: PetscCall(PetscBagGetData(user->param, ¶m));
590: PetscScalar constants[1];
592: constants[0] = param->k;
594: PetscCall(PetscDSSetConstants(ds, 1, constants));
595: }
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: static PetscErrorCode SetupMixedDiscretization(DM dm, AppCtx *user)
600: {
601: DM cdm = dm;
602: PetscFE fe[2];
603: DMPolytopeType ct;
604: PetscInt dim, cStart;
605: MPI_Comm comm;
607: PetscFunctionBeginUser;
608: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
609: PetscCall(DMGetDimension(dm, &dim));
611: /* Create finite element */
612: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
613: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
614: PetscCall(PetscFECreateByCell(comm, dim, dim, ct, "mixed_flux_", PETSC_DEFAULT, &fe[0]));
615: PetscCall(PetscObjectSetName((PetscObject)fe[0], "mixed flux"));
616: /* NOTE: Set the same quadrature order as the primal problem here or use the command line option. */
618: PetscCall(PetscFECreateByCell(comm, dim, 1, ct, "mixed_potential_", PETSC_DEFAULT, &fe[1]));
619: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
620: PetscCall(PetscObjectSetName((PetscObject)fe[1], "mixed potential"));
622: /* Set discretization and boundary conditions for each mesh */
623: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
624: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
625: PetscCall(DMCreateDS(dm));
626: while (cdm) {
627: PetscCall(DMCopyDisc(dm, cdm));
628: PetscCall(DMGetCoarseDM(cdm, &cdm));
629: }
630: PetscCall(PetscFEDestroy(&fe[0]));
631: PetscCall(PetscFEDestroy(&fe[1]));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: PetscErrorCode SetupMixed(DMAdaptor adaptor, DM mdm)
636: {
637: AppCtx *ctx;
639: PetscFunctionBeginUser;
640: PetscCall(DMGetApplicationContext(mdm, &ctx));
641: PetscCall(SetupMixedDiscretization(mdm, ctx));
642: PetscCall(SetupMixedProblem(mdm, ctx));
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: int main(int argc, char **argv)
647: {
648: DM dm, mdm; /* problem specification */
649: SNES snes, msnes; /* nonlinear solvers */
650: Vec u, mu; /* solution vectors */
651: Vec fluxError, exactError; /* Element wise error vector */
652: PetscReal fluxNorm, exactNorm; /* Flux error norm */
653: AppCtx user; /* user-defined work context */
655: PetscFunctionBeginUser;
656: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
657: PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.param));
658: PetscCall(SetupParameters(user.param, &user));
659: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
661: // Set up and solve primal problem
662: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
663: PetscCall(DMSetApplicationContext(dm, &user));
664: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
665: PetscCall(SNESSetDM(snes, dm));
667: PetscCall(SetupPrimalDiscretization(dm, &user));
668: PetscCall(SetupPrimalProblem(dm, &user));
669: PetscCall(DMCreateGlobalVector(dm, &u));
670: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
671: PetscCall(SNESSetFromOptions(snes));
672: PetscCall(DMSNESCheckFromOptions(snes, u));
674: for (PetscInt a = 0; a < user.numAdapt; ++a) {
675: if (a > 0) {
676: char prefix[16];
678: PetscCall(PetscSNPrintf(prefix, 16, "a%d_", (int)a));
679: PetscCall(SNESSetOptionsPrefix(snes, prefix));
680: }
681: PetscCall(SNESSolve(snes, NULL, u));
683: // Needed if you allow SNES to refine
684: PetscCall(SNESGetSolution(snes, &u));
685: PetscCall(VecGetDM(u, &dm));
686: }
688: PetscCall(PetscObjectSetName((PetscObject)u, "Primal Solution "));
689: PetscCall(VecViewFromOptions(u, NULL, "-primal_sol_vec_view"));
691: if (user.byHand) {
692: // Set up and solve mixed problem
693: PetscCall(DMClone(dm, &mdm));
694: PetscCall(SNESCreate(PETSC_COMM_WORLD, &msnes));
695: PetscCall(SNESSetOptionsPrefix(msnes, "mixed_"));
696: PetscCall(SNESSetDM(msnes, mdm));
698: PetscCall(SetupMixedDiscretization(mdm, &user));
699: PetscCall(SetupMixedProblem(mdm, &user));
700: PetscCall(DMCreateGlobalVector(mdm, &mu));
701: PetscCall(DMPlexSetSNESLocalFEM(mdm, PETSC_FALSE, &user));
702: PetscCall(SNESSetFromOptions(msnes));
704: PetscCall(DMSNESCheckFromOptions(msnes, mu));
705: PetscCall(SNESSolve(msnes, NULL, mu));
706: PetscCall(PetscObjectSetName((PetscObject)mu, "Mixed Solution "));
707: PetscCall(VecViewFromOptions(mu, NULL, "-mixed_sol_vec_view"));
709: // Create the error space of piecewise constants
710: DM dmErr;
711: PetscFE feErr;
712: DMPolytopeType ct;
713: PetscInt dim, cStart;
715: PetscCall(DMClone(dm, &dmErr));
716: PetscCall(DMGetDimension(dmErr, &dim));
717: PetscCall(DMPlexGetHeightStratum(dmErr, 0, &cStart, NULL));
718: PetscCall(DMPlexGetCellType(dmErr, cStart, &ct));
719: PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, ct, 0, PETSC_DEFAULT, &feErr));
720: PetscCall(PetscObjectSetName((PetscObject)feErr, "Error"));
721: PetscCall(DMSetField(dmErr, 0, NULL, (PetscObject)feErr));
722: PetscCall(PetscFEDestroy(&feErr));
723: PetscCall(DMCreateDS(dmErr));
724: PetscCall(DMViewFromOptions(dmErr, NULL, "-dmerr_view"));
726: // Compute the flux norm
727: PetscCall(DMGetGlobalVector(dmErr, &fluxError));
728: PetscCall(PetscObjectSetName((PetscObject)fluxError, "Flux Error"));
729: PetscCall(DMGetGlobalVector(dmErr, &exactError));
730: PetscCall(PetscObjectSetName((PetscObject)exactError, "Analytical Error"));
731: PetscCall(DMPlexComputeL2FluxDiffVec(u, 0, mu, 0, fluxError));
732: {
733: PetscDS ds;
734: PetscSimplePointFn *func[2] = {NULL, NULL};
735: void *ctx[2] = {NULL, NULL};
737: PetscCall(DMGetDS(mdm, &ds));
738: PetscCall(PetscDSGetExactSolution(ds, 0, &func[0], &ctx[0]));
739: PetscCall(DMPlexComputeL2DiffVec(mdm, 0.0, func, ctx, mu, exactError));
740: }
741: PetscCall(VecNorm(fluxError, NORM_2, &fluxNorm));
742: PetscCall(VecNorm(exactError, NORM_2, &exactNorm));
743: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Flux error norm = %g\t Exact flux error norm = %g\n", (double)fluxNorm, (double)exactNorm));
744: PetscCall(VecViewFromOptions(fluxError, NULL, "-flux_error_vec_view"));
745: PetscCall(VecViewFromOptions(exactError, NULL, "-exact_error_vec_view"));
747: // Adaptive refinement based on calculated error
748: DM rdm;
749: VecTagger refineTag;
750: DMLabel adaptLabel;
751: IS refineIS;
752: Vec ref;
754: PetscCall(DMLabelCreate(PETSC_COMM_WORLD, "adapt", &adaptLabel));
755: PetscCall(DMLabelSetDefaultValue(adaptLabel, DM_ADAPT_COARSEN));
756: PetscCall(VecTaggerCreate(PETSC_COMM_WORLD, &refineTag));
757: PetscCall(VecTaggerSetFromOptions(refineTag));
758: PetscCall(VecTaggerSetUp(refineTag));
759: PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));
761: PetscCall(VecTaggerComputeIS(refineTag, fluxError, &refineIS, NULL));
762: PetscCall(VecTaggerDestroy(&refineTag));
763: PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
764: PetscCall(ISViewFromOptions(refineIS, NULL, "-refine_is_view"));
765: PetscCall(ISDestroy(&refineIS));
767: PetscCall(DMPlexCreateLabelField(dm, adaptLabel, &ref));
768: PetscCall(VecViewFromOptions(ref, NULL, "-refine_vec_view"));
769: PetscCall(VecDestroy(&ref));
771: // Mark adaptation phase with prefix ref_
772: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, "adapt_"));
773: PetscCall(DMAdaptLabel(dm, adaptLabel, &rdm));
774: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, NULL));
775: PetscCall(PetscObjectSetName((PetscObject)rdm, "Adaptively Refined DM"));
776: PetscCall(DMViewFromOptions(rdm, NULL, "-adapt_dm_view"));
777: PetscCall(DMDestroy(&rdm));
778: PetscCall(DMLabelDestroy(&adaptLabel));
780: // Destroy the error structures
781: PetscCall(DMRestoreGlobalVector(dmErr, &fluxError));
782: PetscCall(DMRestoreGlobalVector(dmErr, &exactError));
783: PetscCall(DMDestroy(&dmErr));
785: // Destroy the mixed structures
786: PetscCall(VecDestroy(&mu));
787: PetscCall(DMDestroy(&mdm));
788: PetscCall(SNESDestroy(&msnes));
789: }
791: // Destroy the primal structures
792: PetscCall(VecDestroy(&u));
793: PetscCall(DMDestroy(&dm));
794: PetscCall(SNESDestroy(&snes));
795: PetscCall(PetscBagDestroy(&user.param));
796: PetscCall(PetscFinalize());
797: return 0;
798: }
800: /*TEST
802: # Tests using the explicit code above
803: testset:
804: suffix: 2d_p2_rt0p0_byhand
805: requires: triangle
806: args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \
807: -primal_potential_petscspace_degree 2 \
808: -mixed_potential_petscdualspace_lagrange_continuity 0 \
809: -mixed_flux_petscspace_type ptrimmed \
810: -mixed_flux_petscspace_components 2 \
811: -mixed_flux_petscspace_ptrimmed_form_degree -1 \
812: -mixed_flux_petscdualspace_order 1 \
813: -mixed_flux_petscdualspace_form_degree -1 \
814: -mixed_flux_petscdualspace_lagrange_trimmed true \
815: -mixed_flux_petscfe_default_quadrature_order 2 \
816: -vec_tagger_type cdf -vec_tagger_box 0.9,1.0 \
817: -tag_view \
818: -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \
819: -dmsnes_check 0.001 -mixed_dmsnes_check 0.001 -pc_type jacobi -mixed_pc_type jacobi
820: test:
821: suffix: quadratic
822: args: -sol_type quadratic
823: test:
824: suffix: trig
825: args: -sol_type trig
826: test:
827: suffix: sensor
828: args: -sol_type sensor
830: # Tests using the embedded adaptor in SNES
831: testset:
832: suffix: 2d_p2_rt0p0
833: requires: triangle defined(PETSC_HAVE_EXECUTABLE_EXPORT)
834: args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \
835: -primal_potential_petscspace_degree 2 \
836: -mixed_potential_petscdualspace_lagrange_continuity 0 \
837: -mixed_flux_petscspace_type ptrimmed \
838: -mixed_flux_petscspace_components 2 \
839: -mixed_flux_petscspace_ptrimmed_form_degree -1 \
840: -mixed_flux_petscdualspace_order 1 \
841: -mixed_flux_petscdualspace_form_degree -1 \
842: -mixed_flux_petscdualspace_lagrange_trimmed true \
843: -mixed_flux_petscfe_default_quadrature_order 2 \
844: -by_hand 0 \
845: -refine_vec_tagger_type cdf -refine_vec_tagger_box 0.9,1.0 \
846: -snes_adapt_view \
847: -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \
848: -adaptor_criterion label -adaptor_type flux -adaptor_mixed_setup_function SetupMixed \
849: -snes_adapt_sequence 1 -pc_type jacobi -mixed_pc_type jacobi
850: test:
851: suffix: quadratic
852: args: -sol_type quadratic -adaptor_monitor_error
853: test:
854: suffix: trig
855: args: -sol_type trig -adaptor_monitor_error
856: test:
857: suffix: sensor
858: args: -sol_type sensor
860: # Tests using multiple adaptor loops
861: testset:
862: suffix: 2d_p2_rt0p0_a2
863: requires: triangle defined(PETSC_HAVE_EXECUTABLE_EXPORT)
864: args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \
865: -primal_potential_petscspace_degree 2 \
866: -mixed_potential_petscdualspace_lagrange_continuity 0 \
867: -mixed_flux_petscspace_type ptrimmed \
868: -mixed_flux_petscspace_components 2 \
869: -mixed_flux_petscspace_ptrimmed_form_degree -1 \
870: -mixed_flux_petscdualspace_order 1 \
871: -mixed_flux_petscdualspace_form_degree -1 \
872: -mixed_flux_petscdualspace_lagrange_trimmed true \
873: -mixed_flux_petscfe_default_quadrature_order 2 \
874: -by_hand 0 \
875: -num_adapt 2 \
876: -refine_vec_tagger_type cdf -refine_vec_tagger_box 0.9,1.0 \
877: -snes_adapt_view \
878: -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \
879: -adaptor_criterion label -adaptor_type gradient -adaptor_mixed_setup_function SetupMixed \
880: -snes_adapt_sequence 2 -pc_type jacobi \
881: -a1_refine_vec_tagger_type cdf -a1_refine_vec_tagger_box 0.9,1.0 \
882: -a1_snes_adapt_view \
883: -a1_adaptor_criterion label -a1_adaptor_type flux -a1_adaptor_mixed_setup_function SetupMixed \
884: -a1_snes_adapt_sequence 1 -a1_pc_type jacobi -a1_mixed_pc_type jacobi
885: test:
886: suffix: sensor
887: args: -sol_type sensor
889: TEST*/