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"; 5 6 /* 7 The primal (or original) Poisson problem, in the strong form, is given by, 8 9 \begin{align} 10 - \nabla \cdot ( \nabla u ) = f 11 \end{align} 12 where $u$ is potential. 13 14 The weak form of this equation is 15 16 \begin{align} 17 < \nabla v, \nabla u > - < v, \nabla u \cdot \hat{n} >_\Gamma - < v, f > = 0 18 \end{align} 19 20 The mixed Poisson problem, in the strong form, is given by, 21 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. 27 28 The weak form of this equation is 29 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} 34 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. 36 37 The following boundary conditions are prescribed. 38 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} 44 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_____________| 59 60 To visualize the automated adaptation 61 62 -dm_adapt_pre_view draw -dm_adapt_view draw -draw_pause -1 -geometry 0,0,1024,1024 63 64 and to compare with a naice gradient estimator use 65 66 -adaptor_type gradient 67 68 To see a sequence of adaptations 69 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 72 73 To get a better view of the by-hand process, use 74 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 81 82 This is also possible with the automated path 83 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 */ 89 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> 97 98 PETSC_EXTERN PetscErrorCode SetupMixed(DMAdaptor, DM); 99 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}; 108 109 typedef struct { 110 PetscBag param; 111 SolType solType; 112 PetscBool byHand; 113 PetscInt numAdapt; 114 } AppCtx; 115 116 typedef struct { 117 PetscReal k; 118 } Parameter; 119 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; 124 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 } 136 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; 141 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; 149 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 } 153 154 /* 155 Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation: 156 157 f:[-1, 1]x[-1, 1] \to R, 158 f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy) 159 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; 167 168 u[0] = PetscSinReal(50. * xy); 169 if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) u[0] *= 0.01; 170 171 return PETSC_SUCCESS; 172 } 173 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; 180 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 } 189 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 } 195 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; 199 200 f0[0] = -2.0 * PetscSqr(n * PETSC_PI) * PetscSinReal(n * PETSC_PI * x[0]) * PetscSinReal(n * PETSC_PI * x[1]); 201 } 202 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; 208 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 } 212 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]); 216 217 for (PetscInt d = 0; d < dim; ++d) f1[d] = k * u_x[uOff_x[0] + d]; 218 } 219 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; 224 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 } 228 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; 233 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 } 237 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]; 242 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 } 246 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]); 250 251 for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = k; 252 } 253 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; 263 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; 272 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 } 277 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 } 282 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]); 286 287 for (PetscInt d = 0; d < dim; d++) f1[d * dim + d] = k * u[uOff[1]]; 288 } 289 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; 294 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 } 298 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; 303 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 } 307 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; 312 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 } 316 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 } 322 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 } 328 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]); 333 334 for (PetscInt d = 0; d < dim; d++) g2[d * dim + d] = k; 335 } 336 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; 344 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 } 351 352 static PetscErrorCode SetupParameters(PetscBag bag, AppCtx *user) 353 { 354 Parameter *param; 355 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 } 363 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 } 374 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; 382 383 PetscFunctionBeginUser; 384 PetscCall(DMGetDS(dm, &ds)); 385 PetscCall(DMGetLabel(dm, "marker", &label)); 386 387 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_primal_uu)); 388 389 switch (user->solType) { 390 case SOL_QUADRATIC: 391 PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_primal, f1_primal)); 392 393 id = 1; 394 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_u, NULL, user, NULL)); 395 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)); 400 401 id = 3; 402 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_u, NULL, user, NULL)); 403 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)); 408 409 PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_u, user)); 410 break; 411 case SOL_TRIG: 412 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_primal, f1_primal)); 413 414 id = 1; 415 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL)); 416 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)); 421 422 id = 3; 423 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL)); 424 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)); 429 430 PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user)); 431 break; 432 case SOL_SENSOR: 433 PetscCall(PetscDSSetResidual(ds, 0, f0_sensor_primal, f1_primal)); 434 435 id = 1; 436 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)sensor_u, NULL, user, NULL)); 437 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)); 442 443 id = 3; 444 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)sensor_u, NULL, user, NULL)); 445 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)); 450 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 } 456 457 /* Setup constants */ 458 { 459 PetscCall(PetscBagGetData(user->param, (void **)¶m)); 460 PetscScalar constants[1]; 461 462 constants[0] = param->k; 463 464 PetscCall(PetscDSSetConstants(ds, 1, constants)); 465 } 466 PetscFunctionReturn(PETSC_SUCCESS); 467 } 468 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; 476 477 PetscFunctionBeginUser; 478 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 479 PetscCall(DMGetDimension(dm, &dim)); 480 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")); 486 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 } 494 495 PetscCall(PetscFEDestroy(&fe[0])); 496 PetscFunctionReturn(PETSC_SUCCESS); 497 } 498 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; 506 507 PetscFunctionBeginUser; 508 PetscCall(DMGetDS(dm, &ds)); 509 PetscCall(DMGetLabel(dm, "marker", &label)); 510 511 /* Residual terms */ 512 PetscCall(PetscDSSetResidual(ds, 0, f0_mixed_q, f1_mixed_q)); 513 514 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_mixed_qq, NULL, NULL, NULL)); 515 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_mixed_qu, NULL)); 516 517 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_mixed_uq, NULL, NULL)); 518 519 switch (user->solType) { 520 case SOL_QUADRATIC: 521 PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_quadratic_u, NULL)); 522 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)); 527 528 id = 2; 529 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_q, NULL, user, NULL)); 530 531 id = 4; 532 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)quadratic_q, NULL, user, NULL)); 533 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)); 538 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)); 544 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)); 549 550 id = 2; 551 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)trig_q, NULL, user, NULL)); 552 553 id = 4; 554 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)trig_q, NULL, user, NULL)); 555 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)); 560 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)); 566 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)); 571 572 id = 2; 573 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)sensor_q, NULL, user, NULL)); 574 575 id = 4; 576 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (PetscVoidFn *)sensor_q, NULL, user, NULL)); 577 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)); 582 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 } 589 590 /* Setup constants */ 591 { 592 PetscCall(PetscBagGetData(user->param, (void **)¶m)); 593 PetscScalar constants[1]; 594 595 constants[0] = param->k; 596 597 PetscCall(PetscDSSetConstants(ds, 1, constants)); 598 } 599 PetscFunctionReturn(PETSC_SUCCESS); 600 } 601 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; 609 610 PetscFunctionBeginUser; 611 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 612 PetscCall(DMGetDimension(dm, &dim)); 613 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. */ 620 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")); 624 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 } 637 638 PetscErrorCode SetupMixed(DMAdaptor adaptor, DM mdm) 639 { 640 AppCtx *ctx; 641 642 PetscFunctionBeginUser; 643 PetscCall(DMGetApplicationContext(mdm, (void **)&ctx)); 644 PetscCall(SetupMixedDiscretization(mdm, ctx)); 645 PetscCall(SetupMixedProblem(mdm, ctx)); 646 PetscFunctionReturn(PETSC_SUCCESS); 647 } 648 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 */ 657 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)); 663 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)); 669 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)); 677 678 for (PetscInt a = 0; a < user.numAdapt; ++a) { 679 if (a > 0) { 680 char prefix[16]; 681 682 PetscCall(PetscSNPrintf(prefix, 16, "a%d_", (int)a)); 683 PetscCall(SNESSetOptionsPrefix(snes, prefix)); 684 } 685 PetscCall(SNESSolve(snes, NULL, u)); 686 687 // Needed if you allow SNES to refine 688 PetscCall(SNESGetSolution(snes, &u)); 689 PetscCall(VecGetDM(u, &dm)); 690 } 691 692 PetscCall(PetscObjectSetName((PetscObject)u, "Primal Solution ")); 693 PetscCall(VecViewFromOptions(u, NULL, "-primal_sol_vec_view")); 694 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)); 701 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)); 708 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")); 713 714 // Create the error space of piecewise constants 715 DM dmErr; 716 PetscFE feErr; 717 DMPolytopeType ct; 718 PetscInt dim, cStart; 719 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")); 730 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}; 741 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")); 751 752 // Adaptive refinement based on calculated error 753 DM rdm; 754 VecTagger refineTag; 755 DMLabel adaptLabel; 756 IS refineIS; 757 Vec ref; 758 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")); 765 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)); 771 772 PetscCall(DMPlexCreateLabelField(dm, adaptLabel, &ref)); 773 PetscCall(VecViewFromOptions(ref, NULL, "-refine_vec_view")); 774 PetscCall(VecDestroy(&ref)); 775 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)); 784 785 // Destroy the error structures 786 PetscCall(DMRestoreGlobalVector(dmErr, &fluxError)); 787 PetscCall(DMRestoreGlobalVector(dmErr, &exactError)); 788 PetscCall(DMDestroy(&dmErr)); 789 790 // Destroy the mixed structures 791 PetscCall(VecDestroy(&mu)); 792 PetscCall(DMDestroy(&mdm)); 793 PetscCall(SNESDestroy(&msnes)); 794 } 795 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 } 804 805 /*TEST 806 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 834 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 864 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 893 894 TEST*/ 895