Actual source code: ex76.c

      1: static char help[] = "Time-dependent Low Mach Flow in 2d and 3d channels with finite elements.\n\
      2: We solve the Low Mach flow problem for both conducting and non-conducting fluids,\n\
      3: using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";

      5: /*F
      6: The non-conducting Low Mach flow is time-dependent isoviscous Navier-Stokes flow. We discretize using the
      7: finite element method on an unstructured mesh. The weak form equations are

      9: \begin{align*}
     10:     < q, \nabla\cdot u > = 0
     11:     <v, du/dt> + <v, u \cdot \nabla u> + < \nabla v, \nu (\nabla u + {\nabla u}^T) > - < \nabla\cdot v, p >  - < v, f  >  = 0
     12:     < w, u \cdot \nabla T > + < \nabla w, \alpha \nabla T > - < w, Q > = 0
     13: \end{align*}

     15: where $\nu$ is the kinematic viscosity and $\alpha$ is thermal diffusivity.

     17: The conducting form is given in the ABLATE documentation [1,2] and derived in Principe and Codina [2].

     19: For visualization, use

     21:   -dm_view hdf5:$PWD/sol.h5 -sol_vec_view hdf5:$PWD/sol.h5::append -exact_vec_view hdf5:$PWD/sol.h5::append

     23: To look at nonlinear solver convergence, use

     25:   -dm_refine <k> -ts_max_steps 1 \
     26:   -ts_view -ts_monitor -snes_monitor -snes_converged_reason -ksp_converged_reason -fieldsplit_pressure_ksp_converged_reason

     28: [1] https://ubchrest.github.io/ablate/content/formulations/lowMachFlow/
     29: [2] https://github.com/UBCHREST/ablate/blob/main/ablateCore/flow/lowMachFlow.c
     30: [3] J. Principe and R. Codina, "Mathematical models for thermally coupled low speed flows", Adv. in Theo. and App. Mech., 2(1), pp.93--112, 2009.
     31: F*/

     33: #include <petscdmplex.h>
     34: #include <petscsnes.h>
     35: #include <petscts.h>
     36: #include <petscds.h>
     37: #include <petscbag.h>

     39: typedef enum {
     40:   MOD_INCOMPRESSIBLE,
     41:   MOD_CONDUCTING,
     42:   NUM_MOD_TYPES
     43: } ModType;
     44: const char *modTypes[NUM_MOD_TYPES + 1] = {"incompressible", "conducting", "unknown"};

     46: typedef enum {
     47:   SOL_QUADRATIC,
     48:   SOL_CUBIC,
     49:   SOL_CUBIC_TRIG,
     50:   SOL_TAYLOR_GREEN,
     51:   SOL_PIPE,
     52:   SOL_PIPE_WIGGLY,
     53:   NUM_SOL_TYPES
     54: } SolType;
     55: const char *solTypes[NUM_SOL_TYPES + 1] = {"quadratic", "cubic", "cubic_trig", "taylor_green", "pipe", "pipe_wiggly", "unknown"};

     57: /* Fields */
     58: const PetscInt VEL  = 0;
     59: const PetscInt PRES = 1;
     60: const PetscInt TEMP = 2;
     61: /* Sources */
     62: const PetscInt MOMENTUM = 0;
     63: const PetscInt MASS     = 1;
     64: const PetscInt ENERGY   = 2;
     65: /* Constants */
     66: const PetscInt STROUHAL = 0;
     67: const PetscInt FROUDE   = 1;
     68: const PetscInt REYNOLDS = 2;
     69: const PetscInt PECLET   = 3;
     70: const PetscInt P_TH     = 4;
     71: const PetscInt MU       = 5;
     72: const PetscInt NU       = 6;
     73: const PetscInt C_P      = 7;
     74: const PetscInt K        = 8;
     75: const PetscInt ALPHA    = 9;
     76: const PetscInt T_IN     = 10;
     77: const PetscInt G_DIR    = 11;
     78: const PetscInt EPSILON  = 12;

     80: typedef struct {
     81:   PetscReal Strouhal; /* Strouhal number */
     82:   PetscReal Froude;   /* Froude number */
     83:   PetscReal Reynolds; /* Reynolds number */
     84:   PetscReal Peclet;   /* Peclet number */
     85:   PetscReal p_th;     /* Thermodynamic pressure */
     86:   PetscReal mu;       /* Dynamic viscosity */
     87:   PetscReal nu;       /* Kinematic viscosity */
     88:   PetscReal c_p;      /* Specific heat at constant pressure */
     89:   PetscReal k;        /* Thermal conductivity */
     90:   PetscReal alpha;    /* Thermal diffusivity */
     91:   PetscReal T_in;     /* Inlet temperature */
     92:   PetscReal g_dir;    /* Gravity direction */
     93:   PetscReal epsilon;  /* Strength of perturbation */
     94: } Parameter;

     96: typedef struct {
     97:   /* Problem definition */
     98:   PetscBag  bag;          /* Holds problem parameters */
     99:   ModType   modType;      /* Model type */
    100:   SolType   solType;      /* MMS solution type */
    101:   PetscBool hasNullSpace; /* Problem has the constant null space for pressure */
    102:   /* Flow diagnostics */
    103:   DM dmCell; /* A DM with piecewise constant discretization */
    104: } AppCtx;

    106: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    107: {
    108:   PetscInt d;
    109:   for (d = 0; d < Nc; ++d) u[d] = 0.0;
    110:   return PETSC_SUCCESS;
    111: }

    113: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
    114: {
    115:   PetscInt d;
    116:   for (d = 0; d < Nc; ++d) u[d] = 1.0;
    117:   return PETSC_SUCCESS;
    118: }

    120: /*
    121:   CASE: quadratic
    122:   In 2D we use exact solution:

    124:     u = t + x^2 + y^2
    125:     v = t + 2x^2 - 2xy
    126:     p = x + y - 1
    127:     T = t + x + y + 1
    128:     f = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2 -4\nu + 2, t (2x - 2y) + 4xy^2 + 2x^2y - 2y^3 -4\nu + 2>
    129:     Q = 1 + 2t + 3x^2 - 2xy + y^2

    131:   so that

    133:     \nabla \cdot u = 2x - 2x = 0

    135:   f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
    136:     = <1, 1> + <t + x^2 + y^2, t + 2x^2 - 2xy> . <<2x, 4x - 2y>, <2y, -2x>> - \nu <4, 4> + <1, 1>
    137:     = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2, t (2x - 2y) + 2x^2y + 4xy^2 - 2y^3> + <-4 \nu + 2, -4\nu + 2>
    138:     = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2 - 4\nu + 2, t (2x - 2y) + 4xy^2 + 2x^2y - 2y^3 - 4\nu + 2>

    140:   Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
    141:     = 1 + <t + x^2 + y^2, t + 2x^2 - 2xy> . <1, 1> - \alpha 0
    142:     = 1 + 2t + 3x^2 - 2xy + y^2
    143: */

    145: static PetscErrorCode quadratic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    146: {
    147:   u[0] = time + X[0] * X[0] + X[1] * X[1];
    148:   u[1] = time + 2.0 * X[0] * X[0] - 2.0 * X[0] * X[1];
    149:   return PETSC_SUCCESS;
    150: }
    151: static PetscErrorCode quadratic_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    152: {
    153:   u[0] = 1.0;
    154:   u[1] = 1.0;
    155:   return PETSC_SUCCESS;
    156: }

    158: static PetscErrorCode quadratic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
    159: {
    160:   p[0] = X[0] + X[1] - 1.0;
    161:   return PETSC_SUCCESS;
    162: }

    164: static PetscErrorCode quadratic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    165: {
    166:   T[0] = time + X[0] + X[1] + 1.0;
    167:   return PETSC_SUCCESS;
    168: }
    169: static PetscErrorCode quadratic_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    170: {
    171:   T[0] = 1.0;
    172:   return PETSC_SUCCESS;
    173: }

    175: static void f0_quadratic_v(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[])
    176: {
    177:   const PetscReal nu = PetscRealPart(constants[NU]);

    179:   f0[0] -= t * (2 * X[0] + 2 * X[1]) + 2 * X[0] * X[0] * X[0] + 4 * X[0] * X[0] * X[1] - 2 * X[0] * X[1] * X[1] - 4.0 * nu + 2;
    180:   f0[1] -= t * (2 * X[0] - 2 * X[1]) + 4 * X[0] * X[1] * X[1] + 2 * X[0] * X[0] * X[1] - 2 * X[1] * X[1] * X[1] - 4.0 * nu + 2;
    181: }

    183: static void f0_quadratic_w(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[])
    184: {
    185:   f0[0] -= 2 * t + 1 + 3 * X[0] * X[0] - 2 * X[0] * X[1] + X[1] * X[1];
    186: }

    188: /*
    189:   CASE: quadratic
    190:   In 2D we use exact solution:

    192:     u = t + x^2 + y^2
    193:     v = t + 2x^2 - 2xy
    194:     p = x + y - 1
    195:     T = t + x + y + 1
    196:   rho = p^{th} / T

    198:   so that

    200:     \nabla \cdot u = 2x - 2x = 0
    201:     grad u = <<2 x, 4x - 2y>, <2 y, -2x>>
    202:     epsilon(u) = 1/2 (grad u + grad u^T) = <<2x, 2x>, <2x, -2x>>
    203:     epsilon'(u) = epsilon(u) - 1/3 (div u) I = epsilon(u)
    204:     div epsilon'(u) = <2, 2>

    206:   f = rho S du/dt + rho u \cdot \nabla u - 2\mu/Re div \epsilon'(u) + \nabla p + rho / F^2 \hat y
    207:     = rho S <1, 1> + rho <t + x^2 + y^2, t + 2x^2 - 2xy> . <<2x, 4x - 2y>, <2y, -2x>> - 2\mu/Re <2, 2> + <1, 1> + rho/F^2 <0, 1>
    208:     = rho S <1, 1> + rho <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2, t (2x - 2y) + 2x^2y + 4xy^2 - 2y^3> - mu/Re <4, 4> + <1, 1> + rho/F^2 <0, 1>

    210:   g = S rho_t + div (rho u)
    211:     = -S pth T_t/T^2 + rho div (u) + u . grad rho
    212:     = -S pth 1/T^2 - pth u . grad T / T^2
    213:     = -pth / T^2 (S + 2t + 3 x^2 - 2xy + y^2)

    215:   Q = rho c_p S dT/dt + rho c_p u . grad T - 1/Pe div k grad T
    216:     = c_p S pth / T + c_p pth (2t + 3 x^2 - 2xy + y^2) / T - k/Pe 0
    217:     = c_p pth / T (S + 2t + 3 x^2 - 2xy + y^2)
    218: */
    219: static void f0_conduct_quadratic_v(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[])
    220: {
    221:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    222:   const PetscReal F    = PetscRealPart(constants[FROUDE]);
    223:   const PetscReal Re   = PetscRealPart(constants[REYNOLDS]);
    224:   const PetscReal mu   = PetscRealPart(constants[MU]);
    225:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    226:   const PetscReal rho  = p_th / (t + X[0] + X[1] + 1.);
    227:   const PetscInt  gd   = (PetscInt)PetscRealPart(constants[G_DIR]);

    229:   f0[0] -= rho * S + rho * (2. * t * (X[0] + X[1]) + 2. * X[0] * X[0] * X[0] + 4. * X[0] * X[0] * X[1] - 2. * X[0] * X[1] * X[1]) - 4. * mu / Re + 1.;
    230:   f0[1] -= rho * S + rho * (2. * t * (X[0] - X[1]) + 2. * X[0] * X[0] * X[1] + 4. * X[0] * X[1] * X[1] - 2. * X[1] * X[1] * X[1]) - 4. * mu / Re + 1.;
    231:   f0[gd] -= rho / PetscSqr(F);
    232: }

    234: static void f0_conduct_quadratic_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[])
    235: {
    236:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    237:   const PetscReal p_th = PetscRealPart(constants[P_TH]);

    239:   f0[0] += p_th * (S + 2. * t + 3. * X[0] * X[0] - 2. * X[0] * X[1] + X[1] * X[1]) / PetscSqr(t + X[0] + X[1] + 1.);
    240: }

    242: static void f0_conduct_quadratic_w(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[])
    243: {
    244:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    245:   const PetscReal c_p  = PetscRealPart(constants[C_P]);
    246:   const PetscReal p_th = PetscRealPart(constants[P_TH]);

    248:   f0[0] -= c_p * p_th * (S + 2. * t + 3. * X[0] * X[0] - 2. * X[0] * X[1] + X[1] * X[1]) / (t + X[0] + X[1] + 1.);
    249: }

    251: /*
    252:   CASE: cubic
    253:   In 2D we use exact solution:

    255:     u = t + x^3 + y^3
    256:     v = t + 2x^3 - 3x^2y
    257:     p = 3/2 x^2 + 3/2 y^2 - 1
    258:     T = t + 1/2 x^2 + 1/2 y^2
    259:     f = < t(3x^2 + 3y^2) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y) + 3x + 1,
    260:           t(3x^2 - 6xy) + 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y + 1>
    261:     Q = x^4 + xy^3 + 2x^3y - 3x^2y^2 + xt + yt - 2\alpha + 1

    263:   so that

    265:     \nabla \cdot u = 3x^2 - 3x^2 = 0

    267:   du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p - f
    268:   = <1,1> + <t(3x^2 + 3y^2) + 3x^5 + 6x^3y^2 - 6x^2y^3, t(3x^2 - 6xy) + 6x^2y^3 + 3x^4y - 6xy^4> - \nu<6x + 6y, 12x - 6y> + <3x, 3y> - <t(3x^2 + 3y^2) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y) + 3x + 1, t(3x^2 - 6xy) + 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y + 1>  = 0

    270:   dT/dt + u \cdot \nabla T - \alpha \Delta T - Q = 1 + (x^3 + y^3) x + (2x^3 - 3x^2y) y - 2*\alpha - (x^4 + xy^3 + 2x^3y - 3x^2y^2 - 2*\alpha +1)   = 0
    271: */
    272: static PetscErrorCode cubic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    273: {
    274:   u[0] = time + X[0] * X[0] * X[0] + X[1] * X[1] * X[1];
    275:   u[1] = time + 2.0 * X[0] * X[0] * X[0] - 3.0 * X[0] * X[0] * X[1];
    276:   return PETSC_SUCCESS;
    277: }
    278: static PetscErrorCode cubic_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    279: {
    280:   u[0] = 1.0;
    281:   u[1] = 1.0;
    282:   return PETSC_SUCCESS;
    283: }

    285: static PetscErrorCode cubic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
    286: {
    287:   p[0] = 3.0 * X[0] * X[0] / 2.0 + 3.0 * X[1] * X[1] / 2.0 - 1.0;
    288:   return PETSC_SUCCESS;
    289: }

    291: static PetscErrorCode cubic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    292: {
    293:   T[0] = time + X[0] * X[0] / 2.0 + X[1] * X[1] / 2.0;
    294:   return PETSC_SUCCESS;
    295: }
    296: static PetscErrorCode cubic_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    297: {
    298:   T[0] = 1.0;
    299:   return PETSC_SUCCESS;
    300: }

    302: static void f0_cubic_v(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[])
    303: {
    304:   const PetscReal nu = PetscRealPart(constants[NU]);

    306:   f0[0] -= (t * (3 * X[0] * X[0] + 3 * X[1] * X[1]) + 3 * X[0] * X[0] * X[0] * X[0] * X[0] + 6 * X[0] * X[0] * X[0] * X[1] * X[1] - 6 * X[0] * X[0] * X[1] * X[1] * X[1] - (6 * X[0] + 6 * X[1]) * nu + 3 * X[0] + 1);
    307:   f0[1] -= (t * (3 * X[0] * X[0] - 6 * X[0] * X[1]) + 3 * X[0] * X[0] * X[0] * X[0] * X[1] + 6 * X[0] * X[0] * X[1] * X[1] * X[1] - 6 * X[0] * X[1] * X[1] * X[1] * X[1] - (12 * X[0] - 6 * X[1]) * nu + 3 * X[1] + 1);
    308: }

    310: static void f0_cubic_w(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[])
    311: {
    312:   const PetscReal alpha = PetscRealPart(constants[ALPHA]);

    314:   f0[0] -= X[0] * X[0] * X[0] * X[0] + 2.0 * X[0] * X[0] * X[0] * X[1] - 3.0 * X[0] * X[0] * X[1] * X[1] + X[0] * X[1] * X[1] * X[1] + X[0] * t + X[1] * t - 2.0 * alpha + 1;
    315: }

    317: /*
    318:   CASE: cubic-trigonometric
    319:   In 2D we use exact solution:

    321:     u = beta cos t + x^3 + y^3
    322:     v = beta sin t + 2x^3 - 3x^2y
    323:     p = 3/2 x^2 + 3/2 y^2 - 1
    324:     T = 20 cos t + 1/2 x^2 + 1/2 y^2
    325:     f = < beta cos t 3x^2         + beta sin t (3y^2 - 1) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y)  + 3x,
    326:           beta cos t (6x^2 - 6xy) - beta sin t (3x^2)     + 3x^4y + 6x^2y^3 - 6xy^4  - \nu(12x - 6y) + 3y>
    327:     Q = beta cos t x + beta sin t (y - 1) + x^4 + 2x^3y - 3x^2y^2 + xy^3 - 2\alpha

    329:   so that

    331:     \nabla \cdot u = 3x^2 - 3x^2 = 0

    333:   f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
    334:     = <-sin t, cos t> + <cos t + x^3 + y^3, sin t + 2x^3 - 3x^2y> <<3x^2, 6x^2 - 6xy>, <3y^2, -3x^2>> - \nu <6x + 6y, 12x - 6y> + <3x, 3y>
    335:     = <-sin t, cos t> + <cos t 3x^2 + 3x^5 + 3x^2y^3 + sin t 3y^2 + 6x^3y^2 - 9x^2y^3, cos t (6x^2 - 6xy) + 6x^5 - 6x^4y + 6x^2y^3 - 6xy^4 + sin t (-3x^2) - 6x^5 + 9x^4y> - \nu <6x + 6y, 12x - 6y> + <3x, 3y>
    336:     = <cos t (3x^2)       + sin t (3y^2 - 1) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu (6x + 6y)  + 3x,
    337:        cos t (6x^2 - 6xy) - sin t (3x^2)     + 3x^4y + 6x^2y^3 - 6xy^4  - \nu (12x - 6y) + 3y>

    339:   Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
    340:     = -sin t + <cos t + x^3 + y^3, sin t + 2x^3 - 3x^2y> . <x, y> - 2 \alpha
    341:     = -sin t + cos t (x) + x^4 + xy^3 + sin t (y) + 2x^3y - 3x^2y^2 - 2 \alpha
    342:     = cos t x + sin t (y - 1) + (x^4 + 2x^3y - 3x^2y^2 + xy^3 - 2 \alpha)
    343: */
    344: static PetscErrorCode cubic_trig_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    345: {
    346:   u[0] = 100. * PetscCosReal(time) + X[0] * X[0] * X[0] + X[1] * X[1] * X[1];
    347:   u[1] = 100. * PetscSinReal(time) + 2.0 * X[0] * X[0] * X[0] - 3.0 * X[0] * X[0] * X[1];
    348:   return PETSC_SUCCESS;
    349: }
    350: static PetscErrorCode cubic_trig_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    351: {
    352:   u[0] = -100. * PetscSinReal(time);
    353:   u[1] = 100. * PetscCosReal(time);
    354:   return PETSC_SUCCESS;
    355: }

    357: static PetscErrorCode cubic_trig_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
    358: {
    359:   p[0] = 3.0 * X[0] * X[0] / 2.0 + 3.0 * X[1] * X[1] / 2.0 - 1.0;
    360:   return PETSC_SUCCESS;
    361: }

    363: static PetscErrorCode cubic_trig_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    364: {
    365:   T[0] = 100. * PetscCosReal(time) + X[0] * X[0] / 2.0 + X[1] * X[1] / 2.0;
    366:   return PETSC_SUCCESS;
    367: }
    368: static PetscErrorCode cubic_trig_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    369: {
    370:   T[0] = -100. * PetscSinReal(time);
    371:   return PETSC_SUCCESS;
    372: }

    374: static void f0_cubic_trig_v(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[])
    375: {
    376:   const PetscReal nu = PetscRealPart(constants[NU]);

    378:   f0[0] -= 100. * PetscCosReal(t) * (3 * X[0] * X[0]) + 100. * PetscSinReal(t) * (3 * X[1] * X[1] - 1.) + 3 * X[0] * X[0] * X[0] * X[0] * X[0] + 6 * X[0] * X[0] * X[0] * X[1] * X[1] - 6 * X[0] * X[0] * X[1] * X[1] * X[1] - (6 * X[0] + 6 * X[1]) * nu + 3 * X[0];
    379:   f0[1] -= 100. * PetscCosReal(t) * (6 * X[0] * X[0] - 6 * X[0] * X[1]) - 100. * PetscSinReal(t) * (3 * X[0] * X[0]) + 3 * X[0] * X[0] * X[0] * X[0] * X[1] + 6 * X[0] * X[0] * X[1] * X[1] * X[1] - 6 * X[0] * X[1] * X[1] * X[1] * X[1] - (12 * X[0] - 6 * X[1]) * nu + 3 * X[1];
    380: }

    382: static void f0_cubic_trig_w(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[])
    383: {
    384:   const PetscReal alpha = PetscRealPart(constants[ALPHA]);

    386:   f0[0] -= 100. * PetscCosReal(t) * X[0] + 100. * PetscSinReal(t) * (X[1] - 1.) + X[0] * X[0] * X[0] * X[0] + 2.0 * X[0] * X[0] * X[0] * X[1] - 3.0 * X[0] * X[0] * X[1] * X[1] + X[0] * X[1] * X[1] * X[1] - 2.0 * alpha;
    387: }

    389: /*
    390:   CASE: Taylor-Green vortex
    391:   In 2D we use exact solution:

    393:     u = 1 - cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)
    394:     v = 1 + sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)
    395:     p = -1/4 [cos(2 \pi(x - t)) + cos(2 \pi(y - t))] exp(-4 \pi^2 \nu t)
    396:     T = t + x + y
    397:     f = <\nu \pi^2 exp(-2\nu \pi^2 t) cos(\pi(x-t)) sin(\pi(y-t)), -\nu \pi^2 exp(-2\nu \pi^2 t) sin(\pi(x-t)) cos(\pi(y-t))  >
    398:     Q = 3 + sin(\pi(x-y)) exp(-2\nu \pi^2 t)

    400:   so that

    402:   \nabla \cdot u = \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t) - \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t) = 0

    404:   f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
    405:     = <-\pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t)) - 2\pi cos(\pi(x - t)) sin(\pi(y - t))) exp(-2 \pi^2 \nu t),
    406:         \pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t)) - 2\pi sin(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t)>
    407:     + < \pi (1 - cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)) sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
    408:         \pi (1 - cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)) cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
    409:     + <-\pi (1 + sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)) cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t),
    410:        -\pi (1 + sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)) sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)>
    411:     + <-2\pi^2 cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
    412:         2\pi^2 sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
    413:     + < \pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t),
    414:         \pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)>
    415:     = <-\pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t)) - 2\pi cos(\pi(x - t)) sin(\pi(y - t))) exp(-2 \pi^2 \nu t),
    416:         \pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t)) - 2\pi sin(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t)>
    417:     + < \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
    418:         \pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
    419:     + <-\pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t),
    420:        -\pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)>
    421:     + <-\pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t),
    422:        -\pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)>
    423:     + <-2\pi^2 cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
    424:         2\pi^2 sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
    425:     + < \pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t),
    426:         \pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)>
    427:     = <-\pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t),
    428:         \pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t)>
    429:     + < \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
    430:         \pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
    431:     + <-\pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t),
    432:        -\pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)>
    433:     = < \pi cos(\pi(x - t)) cos(\pi(y - t)),
    434:         \pi sin(\pi(x - t)) sin(\pi(y - t))>
    435:     + <-\pi cos(\pi(x - t)) cos(\pi(y - t)),
    436:        -\pi sin(\pi(x - t)) sin(\pi(y - t))> = 0
    437:   Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
    438:     = 1 + u \cdot <1, 1> - 0
    439:     = 1 + u + v
    440: */

    442: static PetscErrorCode taylor_green_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    443: {
    444:   u[0] = 1 - PetscCosReal(PETSC_PI * (X[0] - time)) * PetscSinReal(PETSC_PI * (X[1] - time)) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
    445:   u[1] = 1 + PetscSinReal(PETSC_PI * (X[0] - time)) * PetscCosReal(PETSC_PI * (X[1] - time)) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
    446:   return PETSC_SUCCESS;
    447: }
    448: static PetscErrorCode taylor_green_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    449: {
    450:   u[0] = -PETSC_PI * (PetscSinReal(PETSC_PI * (X[0] - time)) * PetscSinReal(PETSC_PI * (X[1] - time)) - PetscCosReal(PETSC_PI * (X[0] - time)) * PetscCosReal(PETSC_PI * (X[1] - time)) - 2 * PETSC_PI * PetscCosReal(PETSC_PI * (X[0] - time)) * PetscSinReal(PETSC_PI * (X[1] - time))) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
    451:   u[1] = PETSC_PI * (PetscSinReal(PETSC_PI * (X[0] - time)) * PetscSinReal(PETSC_PI * (X[1] - time)) - PetscCosReal(PETSC_PI * (X[0] - time)) * PetscCosReal(PETSC_PI * (X[1] - time)) - 2 * PETSC_PI * PetscSinReal(PETSC_PI * (X[0] - time)) * PetscCosReal(PETSC_PI * (X[1] - time))) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
    452:   return PETSC_SUCCESS;
    453: }

    455: static PetscErrorCode taylor_green_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
    456: {
    457:   p[0] = -0.25 * (PetscCosReal(2 * PETSC_PI * (X[0] - time)) + PetscCosReal(2 * PETSC_PI * (X[1] - time))) * PetscExpReal(-4 * PETSC_PI * PETSC_PI * time);
    458:   return PETSC_SUCCESS;
    459: }

    461: static PetscErrorCode taylor_green_p_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
    462: {
    463:   p[0] = PETSC_PI * (0.5 * (PetscSinReal(2 * PETSC_PI * (X[0] - time)) + PetscSinReal(2 * PETSC_PI * (X[1] - time))) + PETSC_PI * (PetscCosReal(2 * PETSC_PI * (X[0] - time)) + PetscCosReal(2 * PETSC_PI * (X[1] - time)))) * PetscExpReal(-4 * PETSC_PI * PETSC_PI * time);
    464:   return PETSC_SUCCESS;
    465: }

    467: static PetscErrorCode taylor_green_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    468: {
    469:   T[0] = time + X[0] + X[1];
    470:   return PETSC_SUCCESS;
    471: }
    472: static PetscErrorCode taylor_green_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    473: {
    474:   T[0] = 1.0;
    475:   return PETSC_SUCCESS;
    476: }

    478: static void f0_taylor_green_w(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[])
    479: {
    480:   PetscScalar vel[2];

    482:   PetscCallAbort(PETSC_COMM_SELF, taylor_green_u(dim, t, X, Nf, vel, NULL));
    483:   f0[0] -= 1.0 + vel[0] + vel[1];
    484: }

    486: /*
    487:   CASE: Pipe flow
    488:   Poiseuille flow, with the incoming fluid having a parabolic temperature profile and the side walls being held at T_in

    490:     u = \Delta Re/(2 mu) y (1 - y)
    491:     v = 0
    492:     p = -\Delta x
    493:     T = y (1 - y) + T_in
    494:   rho = p^{th} / T

    496:   so that

    498:     \nabla \cdot u = 0 - 0 = 0
    499:     grad u = \Delta Re/(2 mu) <<0, 0>, <1 - 2y, 0>>
    500:     epsilon(u) = 1/2 (grad u + grad u^T) = \Delta Re/(4 mu) <<0, 1 - 2y>, <<1 - 2y, 0>>
    501:     epsilon'(u) = epsilon(u) - 1/3 (div u) I = epsilon(u)
    502:     div epsilon'(u) = -\Delta Re/(2 mu) <1, 0>

    504:   f = rho S du/dt + rho u \cdot \nabla u - 2\mu/Re div \epsilon'(u) + \nabla p + rho / F^2 \hat y
    505:     = 0 + 0 - div (2\mu/Re \epsilon'(u) - pI) + rho / F^2 \hat y
    506:     = -\Delta div <<x, (1 - 2y)/2>, <<(1 - 2y)/2, x>> + rho / F^2 \hat y
    507:     = \Delta <1, 0> - \Delta <1, 0> + rho/F^2 <0, 1>
    508:     = rho/F^2 <0, 1>

    510:   g = S rho_t + div (rho u)
    511:     = 0 + rho div (u) + u . grad rho
    512:     = 0 + 0 - pth u . grad T / T^2
    513:     = 0

    515:   Q = rho c_p S dT/dt + rho c_p u . grad T - 1/Pe div k grad T
    516:     = 0 + c_p pth / T 0 + 2 k/Pe
    517:     = 2 k/Pe

    519:   The boundary conditions on the top and bottom are zero velocity and T_in temperature. The boundary term is

    521:     (2\mu/Re \epsilon'(u) - p I) . n = \Delta <<x, (1 - 2y)/2>, <<(1 - 2y)/2, x>> . n

    523:   so that

    525:     x = 0: \Delta <<0, (1 - 2y)/2>, <<(1 - 2y)/2, 0>> . <-1, 0> = <0, (2y - 1)/2>
    526:     x = 1: \Delta <<1, (1 - 2y)/2>, <<(1 - 2y)/2, 1>> . <1, 0> = <1, (1 - 2y)/2>
    527: */
    528: static PetscErrorCode pipe_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    529: {
    530:   Parameter *param = (Parameter *)ctx;

    532:   u[0] = (0.5 * param->Reynolds / param->mu) * X[1] * (1.0 - X[1]);
    533:   u[1] = 0.0;
    534:   return PETSC_SUCCESS;
    535: }
    536: static PetscErrorCode pipe_u_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    537: {
    538:   u[0] = 0.0;
    539:   u[1] = 0.0;
    540:   return PETSC_SUCCESS;
    541: }

    543: static PetscErrorCode pipe_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
    544: {
    545:   p[0] = -X[0];
    546:   return PETSC_SUCCESS;
    547: }
    548: static PetscErrorCode pipe_p_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
    549: {
    550:   p[0] = 0.0;
    551:   return PETSC_SUCCESS;
    552: }

    554: static PetscErrorCode pipe_T(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    555: {
    556:   Parameter *param = (Parameter *)ctx;

    558:   T[0] = X[1] * (1.0 - X[1]) + param->T_in;
    559:   return PETSC_SUCCESS;
    560: }
    561: static PetscErrorCode pipe_T_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    562: {
    563:   T[0] = 0.0;
    564:   return PETSC_SUCCESS;
    565: }

    567: static void f0_conduct_pipe_v(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[])
    568: {
    569:   const PetscReal F    = PetscRealPart(constants[FROUDE]);
    570:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    571:   const PetscReal T_in = PetscRealPart(constants[T_IN]);
    572:   const PetscReal rho  = p_th / (X[1] * (1. - X[1]) + T_in);
    573:   const PetscInt  gd   = (PetscInt)PetscRealPart(constants[G_DIR]);

    575:   f0[gd] -= rho / PetscSqr(F);
    576: }

    578: static void f0_conduct_bd_pipe_v(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[])
    579: {
    580:   PetscReal sigma[4] = {X[0], (PetscReal)(0.5 * (1. - 2. * X[1])), (PetscReal)(0.5 * (1. - 2. * X[1])), X[0]};
    581:   PetscInt  d, e;

    583:   for (d = 0; d < dim; ++d) {
    584:     for (e = 0; e < dim; ++e) f0[d] -= sigma[d * dim + e] * n[e];
    585:   }
    586: }

    588: static void f0_conduct_pipe_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[])
    589: {
    590:   f0[0] += 0.0;
    591: }

    593: static void f0_conduct_pipe_w(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[])
    594: {
    595:   const PetscReal k  = PetscRealPart(constants[K]);
    596:   const PetscReal Pe = PetscRealPart(constants[PECLET]);

    598:   f0[0] -= 2 * k / Pe;
    599: }

    601: /*
    602:   CASE: Wiggly pipe flow
    603:   Perturbed Poiseuille flow, with the incoming fluid having a perturbed parabolic temperature profile and the side walls being held at T_in

    605:     u = \Delta Re/(2 mu) [y (1 - y) + a sin(pi y)]
    606:     v = 0
    607:     p = -\Delta x
    608:     T = y (1 - y) + a sin(pi y) + T_in
    609:   rho = p^{th} / T

    611:   so that

    613:     \nabla \cdot u = 0 - 0 = 0
    614:     grad u = \Delta Re/(2 mu) <<0, 0>, <1 - 2y + a pi cos(pi y), 0>>
    615:     epsilon(u) = 1/2 (grad u + grad u^T) = \Delta Re/(4 mu) <<0, 1 - 2y + a pi cos(pi y)>, <<1 - 2y + a pi cos(pi y), 0>>
    616:     epsilon'(u) = epsilon(u) - 1/3 (div u) I = epsilon(u)
    617:     div epsilon'(u) = -\Delta Re/(2 mu) <1 + a pi^2/2 sin(pi y), 0>

    619:   f = rho S du/dt + rho u \cdot \nabla u - 2\mu/Re div \epsilon'(u) + \nabla p + rho / F^2 \hat y
    620:     = 0 + 0 - div (2\mu/Re \epsilon'(u) - pI) + rho / F^2 \hat y
    621:     = -\Delta div <<x, (1 - 2y)/2 + a pi/2 cos(pi y)>, <<(1 - 2y)/2 + a pi/2 cos(pi y), x>> + rho / F^2 \hat y
    622:     = -\Delta <1 - 1 - a pi^2/2 sin(pi y), 0> + rho/F^2 <0, 1>
    623:     = a \Delta pi^2/2 sin(pi y) <1, 0> + rho/F^2 <0, 1>

    625:   g = S rho_t + div (rho u)
    626:     = 0 + rho div (u) + u . grad rho
    627:     = 0 + 0 - pth u . grad T / T^2
    628:     = 0

    630:   Q = rho c_p S dT/dt + rho c_p u . grad T - 1/Pe div k grad T
    631:     = 0 + c_p pth / T 0 - k/Pe div <0, 1 - 2y + a pi cos(pi y)>
    632:     = - k/Pe (-2 - a pi^2 sin(pi y))
    633:     = 2 k/Pe (1 + a pi^2/2 sin(pi y))

    635:   The boundary conditions on the top and bottom are zero velocity and T_in temperature. The boundary term is

    637:     (2\mu/Re \epsilon'(u) - p I) . n = \Delta <<x, (1 - 2y)/2 + a pi/2 cos(pi y)>, <<(1 - 2y)/2 + a pi/2 cos(pi y), x>> . n

    639:   so that

    641:     x = 0: \Delta <<0, (1 - 2y)/2>, <<(1 - 2y)/2, 0>> . <-1, 0> = <0, (2y - 1)/2 - a pi/2 cos(pi y)>
    642:     x = 1: \Delta <<1, (1 - 2y)/2>, <<(1 - 2y)/2, 1>> . < 1, 0> = <1, (1 - 2y)/2 + a pi/2 cos(pi y)>
    643: */
    644: static PetscErrorCode pipe_wiggly_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    645: {
    646:   Parameter *param = (Parameter *)ctx;

    648:   u[0] = (0.5 * param->Reynolds / param->mu) * (X[1] * (1.0 - X[1]) + param->epsilon * PetscSinReal(PETSC_PI * X[1]));
    649:   u[1] = 0.0;
    650:   return PETSC_SUCCESS;
    651: }
    652: static PetscErrorCode pipe_wiggly_u_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
    653: {
    654:   u[0] = 0.0;
    655:   u[1] = 0.0;
    656:   return PETSC_SUCCESS;
    657: }

    659: static PetscErrorCode pipe_wiggly_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
    660: {
    661:   p[0] = -X[0];
    662:   return PETSC_SUCCESS;
    663: }
    664: static PetscErrorCode pipe_wiggly_p_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
    665: {
    666:   p[0] = 0.0;
    667:   return PETSC_SUCCESS;
    668: }

    670: static PetscErrorCode pipe_wiggly_T(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    671: {
    672:   Parameter *param = (Parameter *)ctx;

    674:   T[0] = X[1] * (1.0 - X[1]) + param->epsilon * PetscSinReal(PETSC_PI * X[1]) + param->T_in;
    675:   return PETSC_SUCCESS;
    676: }
    677: static PetscErrorCode pipe_wiggly_T_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
    678: {
    679:   T[0] = 0.0;
    680:   return PETSC_SUCCESS;
    681: }

    683: static void f0_conduct_pipe_wiggly_v(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[])
    684: {
    685:   const PetscReal F    = PetscRealPart(constants[FROUDE]);
    686:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    687:   const PetscReal T_in = PetscRealPart(constants[T_IN]);
    688:   const PetscReal eps  = PetscRealPart(constants[EPSILON]);
    689:   const PetscReal rho  = p_th / (X[1] * (1. - X[1]) + T_in);
    690:   const PetscInt  gd   = (PetscInt)PetscRealPart(constants[G_DIR]);

    692:   f0[0] -= eps * 0.5 * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * X[1]);
    693:   f0[gd] -= rho / PetscSqr(F);
    694: }

    696: static void f0_conduct_bd_pipe_wiggly_v(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[])
    697: {
    698:   const PetscReal eps      = PetscRealPart(constants[EPSILON]);
    699:   PetscReal       sigma[4] = {X[0], (PetscReal)(0.5 * (1. - 2. * X[1]) + eps * 0.5 * PETSC_PI * PetscCosReal(PETSC_PI * X[1])), (PetscReal)(0.5 * (1. - 2. * X[1]) + eps * 0.5 * PETSC_PI * PetscCosReal(PETSC_PI * X[1])), X[0]};
    700:   PetscInt        d, e;

    702:   for (d = 0; d < dim; ++d) {
    703:     for (e = 0; e < dim; ++e) f0[d] -= sigma[d * dim + e] * n[e];
    704:   }
    705: }

    707: static void f0_conduct_pipe_wiggly_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[])
    708: {
    709:   f0[0] += 0.0;
    710: }

    712: static void f0_conduct_pipe_wiggly_w(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[])
    713: {
    714:   const PetscReal k   = PetscRealPart(constants[K]);
    715:   const PetscReal Pe  = PetscRealPart(constants[PECLET]);
    716:   const PetscReal eps = PetscRealPart(constants[EPSILON]);

    718:   f0[0] -= 2 * k / Pe * (1.0 + eps * 0.5 * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * X[1]));
    719: }

    721: /*      Physics Kernels      */

    723: static void f0_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[])
    724: {
    725:   PetscInt d;
    726:   for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d * dim + d];
    727: }

    729: /* -\frac{Sp^{th}}{T^2} \frac{\partial T}{\partial t} + \frac{p^{th}}{T} \nabla \cdot \vb{u} - \frac{p^{th}}{T^2} \vb{u} \cdot \nabla T */
    730: static void f0_conduct_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[])
    731: {
    732:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    733:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    734:   PetscInt        d;

    736:   // -\frac{S p^{th}}{T^2} \frac{\partial T}{\partial t}
    737:   f0[0] += -u_t[uOff[TEMP]] * S * p_th / PetscSqr(u[uOff[TEMP]]);

    739:   // \frac{p^{th}}{T} \nabla \cdot \vb{u}
    740:   for (d = 0; d < dim; ++d) f0[0] += p_th / u[uOff[TEMP]] * u_x[uOff_x[VEL] + d * dim + d];

    742:   // - \frac{p^{th}}{T^2} \vb{u} \cdot \nabla T
    743:   for (d = 0; d < dim; ++d) f0[0] -= p_th / (u[uOff[TEMP]] * u[uOff[TEMP]]) * u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d];

    745:   // Add in any fixed source term
    746:   if (NfAux > 0) f0[0] += a[aOff[MASS]];
    747: }

    749: /* \vb{u}_t + \vb{u} \cdot \nabla\vb{u} */
    750: static void f0_v(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[])
    751: {
    752:   const PetscInt Nc = dim;
    753:   PetscInt       c, d;

    755:   for (c = 0; c < Nc; ++c) {
    756:     /* \vb{u}_t */
    757:     f0[c] += u_t[uOff[VEL] + c];
    758:     /* \vb{u} \cdot \nabla\vb{u} */
    759:     for (d = 0; d < dim; ++d) f0[c] += u[uOff[VEL] + d] * u_x[uOff_x[VEL] + c * dim + d];
    760:   }
    761: }

    763: /* \rho S \frac{\partial \vb{u}}{\partial t} + \rho \vb{u} \cdot \nabla \vb{u} + \rho \frac{\hat{\vb{z}}}{F^2} */
    764: static void f0_conduct_v(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[])
    765: {
    766:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    767:   const PetscReal F    = PetscRealPart(constants[FROUDE]);
    768:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    769:   const PetscReal rho  = p_th / PetscRealPart(u[uOff[TEMP]]);
    770:   const PetscInt  gdir = (PetscInt)PetscRealPart(constants[G_DIR]);
    771:   PetscInt        Nc   = dim;
    772:   PetscInt        c, d;

    774:   // \rho S \frac{\partial \vb{u}}{\partial t}
    775:   for (d = 0; d < dim; ++d) f0[d] = rho * S * u_t[uOff[VEL] + d];

    777:   // \rho \vb{u} \cdot \nabla \vb{u}
    778:   for (c = 0; c < Nc; ++c) {
    779:     for (d = 0; d < dim; ++d) f0[c] += rho * u[uOff[VEL] + d] * u_x[uOff_x[VEL] + c * dim + d];
    780:   }

    782:   // rho \hat{z}/F^2
    783:   f0[gdir] += rho / (F * F);

    785:   // Add in any fixed source term
    786:   if (NfAux > 0) {
    787:     for (d = 0; d < dim; ++d) f0[d] += a[aOff[MOMENTUM] + d];
    788:   }
    789: }

    791: /*f1_v = \nu[grad(u) + grad(u)^T] - pI */
    792: static void f1_v(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[])
    793: {
    794:   const PetscReal nu = PetscRealPart(constants[NU]);
    795:   const PetscInt  Nc = dim;
    796:   PetscInt        c, d;

    798:   for (c = 0; c < Nc; ++c) {
    799:     for (d = 0; d < dim; ++d) f1[c * dim + d] = nu * (u_x[c * dim + d] + u_x[d * dim + c]);
    800:     f1[c * dim + c] -= u[uOff[1]];
    801:   }
    802: }

    804: /* 2 \mu/Re (1/2 (\nabla \vb{u} + \nabla \vb{u}^T) - 1/3 (\nabla \cdot \vb{u}) I) - p I */
    805: static void f1_conduct_v(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[])
    806: {
    807:   const PetscReal Re    = PetscRealPart(constants[REYNOLDS]);
    808:   const PetscReal mu    = PetscRealPart(constants[MU]);
    809:   const PetscReal coef  = mu / Re;
    810:   PetscReal       u_div = 0.0;
    811:   const PetscInt  Nc    = dim;
    812:   PetscInt        c, d;

    814:   for (c = 0; c < Nc; ++c) u_div += PetscRealPart(u_x[uOff_x[VEL] + c * dim + c]);

    816:   for (c = 0; c < Nc; ++c) {
    817:     // 2 \mu/Re 1/2 (\nabla \vb{u} + \nabla \vb{u}^T
    818:     for (d = 0; d < dim; ++d) f1[c * dim + d] += coef * (u_x[uOff_x[VEL] + c * dim + d] + u_x[uOff_x[VEL] + d * dim + c]);
    819:     // -2/3 \mu/Re (\nabla \cdot \vb{u}) I
    820:     f1[c * dim + c] -= 2.0 * coef / 3.0 * u_div;
    821:   }

    823:   // -p I
    824:   for (c = 0; c < Nc; ++c) f1[c * dim + c] -= u[uOff[PRES]];
    825: }

    827: /* T_t + \vb{u} \cdot \nabla T */
    828: static void f0_w(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[])
    829: {
    830:   PetscInt d;

    832:   /* T_t */
    833:   f0[0] += u_t[uOff[TEMP]];
    834:   /* \vb{u} \cdot \nabla T */
    835:   for (d = 0; d < dim; ++d) f0[0] += u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d];
    836: }

    838: /* \frac{C_p S p^{th}}{T} \frac{\partial T}{\partial t} + \frac{C_p p^{th}}{T} \vb{u} \cdot \nabla T */
    839: static void f0_conduct_w(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[])
    840: {
    841:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    842:   const PetscReal c_p  = PetscRealPart(constants[C_P]);
    843:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    844:   PetscInt        d;

    846:   // \frac{C_p S p^{th}}{T} \frac{\partial T}{\partial t}
    847:   f0[0] = c_p * S * p_th / u[uOff[TEMP]] * u_t[uOff[TEMP]];

    849:   // \frac{C_p p^{th}}{T} \vb{u} \cdot \nabla T
    850:   for (d = 0; d < dim; ++d) f0[0] += c_p * p_th / u[uOff[TEMP]] * u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d];

    852:   // Add in any fixed source term
    853:   if (NfAux > 0) f0[0] += a[aOff[ENERGY]];
    854: }

    856: static void f1_w(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[])
    857: {
    858:   const PetscReal alpha = PetscRealPart(constants[ALPHA]);
    859:   PetscInt        d;

    861:   for (d = 0; d < dim; ++d) f1[d] = alpha * u_x[uOff_x[2] + d];
    862: }

    864: /* \frac{k}{Pe} \nabla T */
    865: static void f1_conduct_w(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[])
    866: {
    867:   const PetscReal Pe = PetscRealPart(constants[PECLET]);
    868:   const PetscReal k  = PetscRealPart(constants[K]);
    869:   PetscInt        d;

    871:   // \frac{k}{Pe} \nabla T
    872:   for (d = 0; d < dim; ++d) f1[d] = k / Pe * u_x[uOff_x[TEMP] + d];
    873: }

    875: static void g1_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 g1[])
    876: {
    877:   PetscInt d;
    878:   for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
    879: }

    881: static void g0_vu(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[])
    882: {
    883:   PetscInt       c, d;
    884:   const PetscInt Nc = dim;

    886:   for (d = 0; d < dim; ++d) g0[d * dim + d] = u_tShift;

    888:   for (c = 0; c < Nc; ++c) {
    889:     for (d = 0; d < dim; ++d) g0[c * Nc + d] += u_x[c * Nc + d];
    890:   }
    891: }

    893: static void g1_vu(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[])
    894: {
    895:   PetscInt NcI = dim;
    896:   PetscInt NcJ = dim;
    897:   PetscInt c, d, e;

    899:   for (c = 0; c < NcI; ++c) {
    900:     for (d = 0; d < NcJ; ++d) {
    901:       for (e = 0; e < dim; ++e) {
    902:         if (c == d) g1[(c * NcJ + d) * dim + e] += u[e];
    903:       }
    904:     }
    905:   }
    906: }

    908: static void g0_conduct_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 g0[])
    909: {
    910:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    911:   PetscInt        d;

    913:   // - \phi_i \frac{p^{th}}{T^2} \frac{\partial T}{\partial x_c} \psi_{j, u_c}
    914:   for (d = 0; d < dim; ++d) g0[d] = -p_th / PetscSqr(u[uOff[TEMP]]) * u_x[uOff_x[TEMP] + d];
    915: }

    917: static void g1_conduct_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 g1[])
    918: {
    919:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    920:   PetscInt        d;

    922:   // \phi_i \frac{p^{th}}{T} \frac{\partial \psi_{u_c,j}}{\partial x_c}
    923:   for (d = 0; d < dim; ++d) g1[d * dim + d] = p_th / u[uOff[TEMP]];
    924: }

    926: static void g0_conduct_qT(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[])
    927: {
    928:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    929:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    930:   PetscInt        d;

    932:   // - \phi_i \frac{S p^{th}}{T^2} \psi_j
    933:   g0[0] -= S * p_th / PetscSqr(u[uOff[TEMP]]) * u_tShift;
    934:   // \phi_i 2 \frac{S p^{th}}{T^3} T_t \psi_j
    935:   g0[0] += 2.0 * S * p_th / PetscPowScalarInt(u[uOff[TEMP]], 3) * u_t[uOff[TEMP]];
    936:   // \phi_i \frac{p^{th}}{T^2} \left( - \nabla \cdot \vb{u} \psi_j + \frac{2}{T} \vb{u} \cdot \nabla T \psi_j \right)
    937:   for (d = 0; d < dim; ++d) g0[0] += p_th / PetscSqr(u[uOff[TEMP]]) * (-u_x[uOff_x[VEL] + d * dim + d] + 2.0 / u[uOff[TEMP]] * u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d]);
    938: }

    940: static void g1_conduct_qT(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[])
    941: {
    942:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    943:   PetscInt        d;

    945:   // - \phi_i \frac{p^{th}}{T^2} \vb{u} \cdot \nabla \psi_j
    946:   for (d = 0; d < dim; ++d) g1[d] = -p_th / PetscSqr(u[uOff[TEMP]]) * u[uOff[VEL] + d];
    947: }

    949: static void g2_vp(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[])
    950: {
    951:   PetscInt d;
    952:   for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0;
    953: }

    955: static void g3_vu(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[])
    956: {
    957:   const PetscReal nu = PetscRealPart(constants[NU]);
    958:   const PetscInt  Nc = dim;
    959:   PetscInt        c, d;

    961:   for (c = 0; c < Nc; ++c) {
    962:     for (d = 0; d < dim; ++d) {
    963:       g3[((c * Nc + c) * dim + d) * dim + d] += nu;
    964:       g3[((c * Nc + d) * dim + d) * dim + c] += nu;
    965:     }
    966:   }
    967: }

    969: static void g0_conduct_vT(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[])
    970: {
    971:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    972:   const PetscReal F    = PetscRealPart(constants[FROUDE]);
    973:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    974:   const PetscInt  gdir = (PetscInt)PetscRealPart(constants[G_DIR]);
    975:   const PetscInt  Nc   = dim;
    976:   PetscInt        c, d;

    978:   // - \vb{\phi}_i \cdot \vb{u}_t \frac{p^{th} S}{T^2} \psi_j
    979:   for (d = 0; d < dim; ++d) g0[d] -= p_th * S / PetscSqr(u[uOff[TEMP]]) * u_t[uOff[VEL] + d];

    981:   // - \vb{\phi}_i \cdot \vb{u} \cdot \nabla \vb{u} \frac{p^{th}}{T^2} \psi_j
    982:   for (c = 0; c < Nc; ++c) {
    983:     for (d = 0; d < dim; ++d) g0[c] -= p_th / PetscSqr(u[uOff[TEMP]]) * u[uOff[VEL] + d] * u_x[uOff_x[VEL] + c * dim + d];
    984:   }

    986:   // - \vb{\phi}_i \cdot \vu{z} \frac{p^{th}}{T^2 F^2} \psi_j
    987:   g0[gdir] -= p_th / PetscSqr(u[uOff[TEMP]] * F);
    988: }

    990: static void g0_conduct_vu(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[])
    991: {
    992:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    993:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    994:   const PetscInt  Nc   = dim;
    995:   PetscInt        c, d;

    997:   // \vb{\phi}_i \cdot S \rho \psi_j
    998:   for (d = 0; d < dim; ++d) g0[d * dim + d] = S * p_th / u[uOff[TEMP]] * u_tShift;

    1000:   // \phi^c_i \cdot \rho \frac{\partial u^c}{\partial x^d} \psi^d_j
    1001:   for (c = 0; c < Nc; ++c) {
    1002:     for (d = 0; d < dim; ++d) g0[c * Nc + d] += p_th / u[uOff[TEMP]] * u_x[uOff_x[VEL] + c * Nc + d];
    1003:   }
    1004: }

    1006: static void g1_conduct_vu(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[])
    1007: {
    1008:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    1009:   const PetscInt  NcI  = dim;
    1010:   const PetscInt  NcJ  = dim;
    1011:   PetscInt        c, d, e;

    1013:   // \phi^c_i \rho u^e \frac{\partial \psi^d_j}{\partial x^e}
    1014:   for (c = 0; c < NcI; ++c) {
    1015:     for (d = 0; d < NcJ; ++d) {
    1016:       for (e = 0; e < dim; ++e) {
    1017:         if (c == d) g1[(c * NcJ + d) * dim + e] += p_th / u[uOff[TEMP]] * u[uOff[VEL] + e];
    1018:       }
    1019:     }
    1020:   }
    1021: }

    1023: static void g3_conduct_vu(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[])
    1024: {
    1025:   const PetscReal Re = PetscRealPart(constants[REYNOLDS]);
    1026:   const PetscReal mu = PetscRealPart(constants[MU]);
    1027:   const PetscInt  Nc = dim;
    1028:   PetscInt        c, d;

    1030:   for (c = 0; c < Nc; ++c) {
    1031:     for (d = 0; d < dim; ++d) {
    1032:       // \frac{\partial \phi^c_i}{\partial x^d} \mu/Re \frac{\partial \psi^c_i}{\partial x^d}
    1033:       g3[((c * Nc + c) * dim + d) * dim + d] += mu / Re; // gradU
    1034:       // \frac{\partial \phi^c_i}{\partial x^d} \mu/Re \frac{\partial \psi^d_i}{\partial x^c}
    1035:       g3[((c * Nc + d) * dim + d) * dim + c] += mu / Re; // gradU transpose
    1036:       // \frac{\partial \phi^c_i}{\partial x^d} -2/3 \mu/Re \frac{\partial \psi^d_i}{\partial x^c}
    1037:       g3[((c * Nc + d) * dim + c) * dim + d] -= 2.0 / 3.0 * mu / Re;
    1038:     }
    1039:   }
    1040: }

    1042: static void g2_conduct_vp(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[])
    1043: {
    1044:   PetscInt d;
    1045:   for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0;
    1046: }

    1048: static void g0_wT(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[])
    1049: {
    1050:   g0[0] = u_tShift;
    1051: }

    1053: static void g0_wu(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[])
    1054: {
    1055:   PetscInt d;
    1056:   for (d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2] + d];
    1057: }

    1059: static void g1_wT(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[])
    1060: {
    1061:   PetscInt d;
    1062:   for (d = 0; d < dim; ++d) g1[d] = u[uOff[0] + d];
    1063: }

    1065: static void g3_wT(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[])
    1066: {
    1067:   const PetscReal alpha = PetscRealPart(constants[ALPHA]);
    1068:   PetscInt        d;

    1070:   for (d = 0; d < dim; ++d) g3[d * dim + d] = alpha;
    1071: }

    1073: static void g0_conduct_wu(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[])
    1074: {
    1075:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    1076:   const PetscReal c_p  = PetscRealPart(constants[C_P]);
    1077:   PetscInt        d;

    1079:   // \phi_i \frac{C_p p^{th}}{T} \nabla T \cdot \psi_j
    1080:   for (d = 0; d < dim; ++d) g0[d] = c_p * p_th / u[uOff[TEMP]] * u_x[uOff_x[TEMP] + d];
    1081: }

    1083: static void g0_conduct_wT(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[])
    1084: {
    1085:   const PetscReal S    = PetscRealPart(constants[STROUHAL]);
    1086:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    1087:   const PetscReal c_p  = PetscRealPart(constants[C_P]);
    1088:   PetscInt        d;

    1090:   // \psi_i C_p S p^{th}\T \psi_{j}
    1091:   g0[0] += c_p * S * p_th / u[uOff[TEMP]] * u_tShift;
    1092:   // - \phi_i C_p S p^{th}/T^2 T_t \psi_j
    1093:   g0[0] -= c_p * S * p_th / PetscSqr(u[uOff[TEMP]]) * u_t[uOff[TEMP]];
    1094:   // - \phi_i C_p p^{th}/T^2 \vb{u} \cdot \nabla T \psi_j
    1095:   for (d = 0; d < dim; ++d) g0[0] -= c_p * p_th / PetscSqr(u[uOff[TEMP]]) * u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d];
    1096: }

    1098: static void g1_conduct_wT(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[])
    1099: {
    1100:   const PetscReal p_th = PetscRealPart(constants[P_TH]);
    1101:   const PetscReal c_p  = PetscRealPart(constants[C_P]);
    1102:   PetscInt        d;

    1104:   // \phi_i C_p p^{th}/T \vb{u} \cdot \nabla \psi_j
    1105:   for (d = 0; d < dim; ++d) g1[d] += c_p * p_th / u[uOff[TEMP]] * u[uOff[VEL] + d];
    1106: }

    1108: static void g3_conduct_wT(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[])
    1109: {
    1110:   const PetscReal Pe = PetscRealPart(constants[PECLET]);
    1111:   const PetscReal k  = PetscRealPart(constants[K]);
    1112:   PetscInt        d;

    1114:   // \nabla \phi_i \frac{k}{Pe} \nabla \phi_j
    1115:   for (d = 0; d < dim; ++d) g3[d * dim + d] = k / Pe;
    1116: }

    1118: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
    1119: {
    1120:   PetscInt mod, sol;

    1122:   PetscFunctionBeginUser;
    1123:   options->modType      = MOD_INCOMPRESSIBLE;
    1124:   options->solType      = SOL_QUADRATIC;
    1125:   options->hasNullSpace = PETSC_TRUE;
    1126:   options->dmCell       = NULL;

    1128:   PetscOptionsBegin(comm, "", "Low Mach flow Problem Options", "DMPLEX");
    1129:   mod = options->modType;
    1130:   PetscCall(PetscOptionsEList("-mod_type", "The model type", "ex76.c", modTypes, NUM_MOD_TYPES, modTypes[options->modType], &mod, NULL));
    1131:   options->modType = (ModType)mod;
    1132:   sol              = options->solType;
    1133:   PetscCall(PetscOptionsEList("-sol_type", "The solution type", "ex76.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL));
    1134:   options->solType = (SolType)sol;
    1135:   PetscOptionsEnd();
    1136:   PetscFunctionReturn(PETSC_SUCCESS);
    1137: }

    1139: static PetscErrorCode SetupParameters(DM dm, AppCtx *user)
    1140: {
    1141:   PetscBag   bag;
    1142:   Parameter *p;
    1143:   PetscReal  dir;
    1144:   PetscInt   dim;

    1146:   PetscFunctionBeginUser;
    1147:   PetscCall(DMGetDimension(dm, &dim));
    1148:   dir = (PetscReal)(dim - 1);
    1149:   /* setup PETSc parameter bag */
    1150:   PetscCall(PetscBagGetData(user->bag, (void **)&p));
    1151:   PetscCall(PetscBagSetName(user->bag, "par", "Low Mach flow parameters"));
    1152:   bag = user->bag;
    1153:   PetscCall(PetscBagRegisterReal(bag, &p->Strouhal, 1.0, "S", "Strouhal number"));
    1154:   PetscCall(PetscBagRegisterReal(bag, &p->Froude, 1.0, "Fr", "Froude number"));
    1155:   PetscCall(PetscBagRegisterReal(bag, &p->Reynolds, 1.0, "Re", "Reynolds number"));
    1156:   PetscCall(PetscBagRegisterReal(bag, &p->Peclet, 1.0, "Pe", "Peclet number"));
    1157:   PetscCall(PetscBagRegisterReal(bag, &p->p_th, 1.0, "p_th", "Thermodynamic pressure"));
    1158:   PetscCall(PetscBagRegisterReal(bag, &p->mu, 1.0, "mu", "Dynamic viscosity"));
    1159:   PetscCall(PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity"));
    1160:   PetscCall(PetscBagRegisterReal(bag, &p->c_p, 1.0, "c_p", "Specific heat at constant pressure"));
    1161:   PetscCall(PetscBagRegisterReal(bag, &p->k, 1.0, "k", "Thermal conductivity"));
    1162:   PetscCall(PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity"));
    1163:   PetscCall(PetscBagRegisterReal(bag, &p->T_in, 1.0, "T_in", "Inlet temperature"));
    1164:   PetscCall(PetscBagRegisterReal(bag, &p->g_dir, dir, "g_dir", "Gravity direction"));
    1165:   PetscCall(PetscBagRegisterReal(bag, &p->epsilon, 1.0, "epsilon", "Perturbation strength"));
    1166:   PetscFunctionReturn(PETSC_SUCCESS);
    1167: }

    1169: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
    1170: {
    1171:   PetscFunctionBeginUser;
    1172:   PetscCall(DMCreate(comm, dm));
    1173:   PetscCall(DMSetType(*dm, DMPLEX));
    1174:   PetscCall(DMSetFromOptions(*dm));
    1175:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
    1176:   PetscFunctionReturn(PETSC_SUCCESS);
    1177: }

    1179: static PetscErrorCode UniformBoundaryConditions(DM dm, DMLabel label, PetscSimplePointFn *exactFuncs[], PetscSimplePointFn *exactFuncs_t[], AppCtx *user)
    1180: {
    1181:   PetscDS  ds;
    1182:   PetscInt id;
    1183:   void    *ctx;

    1185:   PetscFunctionBeginUser;
    1186:   PetscCall(DMGetDS(dm, &ds));
    1187:   PetscCall(PetscBagGetData(user->bag, &ctx));
    1188:   id = 3;
    1189:   PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
    1190:   id = 1;
    1191:   PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
    1192:   id = 2;
    1193:   PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "right wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
    1194:   id = 4;
    1195:   PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
    1196:   id = 3;
    1197:   PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall temp", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1198:   id = 1;
    1199:   PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall temp", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1200:   id = 2;
    1201:   PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "right wall temp", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1202:   id = 4;
    1203:   PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall temp", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1204:   PetscFunctionReturn(PETSC_SUCCESS);
    1205: }

    1207: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
    1208: {
    1209:   PetscSimplePointFn *exactFuncs[3];
    1210:   PetscSimplePointFn *exactFuncs_t[3];
    1211:   PetscDS             ds;
    1212:   PetscWeakForm       wf;
    1213:   DMLabel             label;
    1214:   Parameter          *ctx;
    1215:   PetscInt            id, bd;

    1217:   PetscFunctionBeginUser;
    1218:   PetscCall(DMGetLabel(dm, "marker", &label));
    1219:   PetscCall(DMGetDS(dm, &ds));
    1220:   PetscCall(PetscDSGetWeakForm(ds, &wf));

    1222:   switch (user->modType) {
    1223:   case MOD_INCOMPRESSIBLE:
    1224:     PetscCall(PetscDSSetResidual(ds, VEL, f0_v, f1_v));
    1225:     PetscCall(PetscDSSetResidual(ds, PRES, f0_q, NULL));
    1226:     PetscCall(PetscDSSetResidual(ds, TEMP, f0_w, f1_w));

    1228:     PetscCall(PetscDSSetJacobian(ds, VEL, VEL, g0_vu, g1_vu, NULL, g3_vu));
    1229:     PetscCall(PetscDSSetJacobian(ds, VEL, PRES, NULL, NULL, g2_vp, NULL));
    1230:     PetscCall(PetscDSSetJacobian(ds, PRES, VEL, NULL, g1_qu, NULL, NULL));
    1231:     PetscCall(PetscDSSetJacobian(ds, TEMP, VEL, g0_wu, NULL, NULL, NULL));
    1232:     PetscCall(PetscDSSetJacobian(ds, TEMP, TEMP, g0_wT, g1_wT, NULL, g3_wT));

    1234:     switch (user->solType) {
    1235:     case SOL_QUADRATIC:
    1236:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_quadratic_v, 0, NULL));
    1237:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_quadratic_w, 0, NULL));

    1239:       exactFuncs[VEL]    = quadratic_u;
    1240:       exactFuncs[PRES]   = quadratic_p;
    1241:       exactFuncs[TEMP]   = quadratic_T;
    1242:       exactFuncs_t[VEL]  = quadratic_u_t;
    1243:       exactFuncs_t[PRES] = NULL;
    1244:       exactFuncs_t[TEMP] = quadratic_T_t;

    1246:       PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
    1247:       break;
    1248:     case SOL_CUBIC:
    1249:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_cubic_v, 0, NULL));
    1250:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_cubic_w, 0, NULL));

    1252:       exactFuncs[VEL]    = cubic_u;
    1253:       exactFuncs[PRES]   = cubic_p;
    1254:       exactFuncs[TEMP]   = cubic_T;
    1255:       exactFuncs_t[VEL]  = cubic_u_t;
    1256:       exactFuncs_t[PRES] = NULL;
    1257:       exactFuncs_t[TEMP] = cubic_T_t;

    1259:       PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
    1260:       break;
    1261:     case SOL_CUBIC_TRIG:
    1262:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_cubic_trig_v, 0, NULL));
    1263:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_cubic_trig_w, 0, NULL));

    1265:       exactFuncs[VEL]    = cubic_trig_u;
    1266:       exactFuncs[PRES]   = cubic_trig_p;
    1267:       exactFuncs[TEMP]   = cubic_trig_T;
    1268:       exactFuncs_t[VEL]  = cubic_trig_u_t;
    1269:       exactFuncs_t[PRES] = NULL;
    1270:       exactFuncs_t[TEMP] = cubic_trig_T_t;

    1272:       PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
    1273:       break;
    1274:     case SOL_TAYLOR_GREEN:
    1275:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_taylor_green_w, 0, NULL));

    1277:       exactFuncs[VEL]    = taylor_green_u;
    1278:       exactFuncs[PRES]   = taylor_green_p;
    1279:       exactFuncs[TEMP]   = taylor_green_T;
    1280:       exactFuncs_t[VEL]  = taylor_green_u_t;
    1281:       exactFuncs_t[PRES] = taylor_green_p_t;
    1282:       exactFuncs_t[TEMP] = taylor_green_T_t;

    1284:       PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
    1285:       break;
    1286:     default:
    1287:       SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
    1288:     }
    1289:     break;
    1290:   case MOD_CONDUCTING:
    1291:     PetscCall(PetscDSSetResidual(ds, VEL, f0_conduct_v, f1_conduct_v));
    1292:     PetscCall(PetscDSSetResidual(ds, PRES, f0_conduct_q, NULL));
    1293:     PetscCall(PetscDSSetResidual(ds, TEMP, f0_conduct_w, f1_conduct_w));

    1295:     PetscCall(PetscDSSetJacobian(ds, VEL, VEL, g0_conduct_vu, g1_conduct_vu, NULL, g3_conduct_vu));
    1296:     PetscCall(PetscDSSetJacobian(ds, VEL, PRES, NULL, NULL, g2_conduct_vp, NULL));
    1297:     PetscCall(PetscDSSetJacobian(ds, VEL, TEMP, g0_conduct_vT, NULL, NULL, NULL));
    1298:     PetscCall(PetscDSSetJacobian(ds, PRES, VEL, g0_conduct_qu, g1_conduct_qu, NULL, NULL));
    1299:     PetscCall(PetscDSSetJacobian(ds, PRES, TEMP, g0_conduct_qT, g1_conduct_qT, NULL, NULL));
    1300:     PetscCall(PetscDSSetJacobian(ds, TEMP, VEL, g0_conduct_wu, NULL, NULL, NULL));
    1301:     PetscCall(PetscDSSetJacobian(ds, TEMP, TEMP, g0_conduct_wT, g1_conduct_wT, NULL, g3_conduct_wT));

    1303:     switch (user->solType) {
    1304:     case SOL_QUADRATIC:
    1305:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_conduct_quadratic_v, 0, NULL));
    1306:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, PRES, 0, 1, f0_conduct_quadratic_q, 0, NULL));
    1307:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_conduct_quadratic_w, 0, NULL));

    1309:       exactFuncs[VEL]    = quadratic_u;
    1310:       exactFuncs[PRES]   = quadratic_p;
    1311:       exactFuncs[TEMP]   = quadratic_T;
    1312:       exactFuncs_t[VEL]  = quadratic_u_t;
    1313:       exactFuncs_t[PRES] = NULL;
    1314:       exactFuncs_t[TEMP] = quadratic_T_t;

    1316:       PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
    1317:       break;
    1318:     case SOL_PIPE:
    1319:       user->hasNullSpace = PETSC_FALSE;
    1320:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_conduct_pipe_v, 0, NULL));
    1321:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, PRES, 0, 1, f0_conduct_pipe_q, 0, NULL));
    1322:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_conduct_pipe_w, 0, NULL));

    1324:       exactFuncs[VEL]    = pipe_u;
    1325:       exactFuncs[PRES]   = pipe_p;
    1326:       exactFuncs[TEMP]   = pipe_T;
    1327:       exactFuncs_t[VEL]  = pipe_u_t;
    1328:       exactFuncs_t[PRES] = pipe_p_t;
    1329:       exactFuncs_t[TEMP] = pipe_T_t;

    1331:       PetscCall(PetscBagGetData(user->bag, (void **)&ctx));
    1332:       id = 2;
    1333:       PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
    1334:       PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    1335:       PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_v, 0, NULL));
    1336:       id = 4;
    1337:       PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
    1338:       PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    1339:       PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_v, 0, NULL));
    1340:       id = 4;
    1341:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1342:       id = 3;
    1343:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
    1344:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1345:       id = 1;
    1346:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
    1347:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1348:       break;
    1349:     case SOL_PIPE_WIGGLY:
    1350:       user->hasNullSpace = PETSC_FALSE;
    1351:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_conduct_pipe_wiggly_v, 0, NULL));
    1352:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, PRES, 0, 1, f0_conduct_pipe_wiggly_q, 0, NULL));
    1353:       PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_conduct_pipe_wiggly_w, 0, NULL));

    1355:       exactFuncs[VEL]    = pipe_wiggly_u;
    1356:       exactFuncs[PRES]   = pipe_wiggly_p;
    1357:       exactFuncs[TEMP]   = pipe_wiggly_T;
    1358:       exactFuncs_t[VEL]  = pipe_wiggly_u_t;
    1359:       exactFuncs_t[PRES] = pipe_wiggly_p_t;
    1360:       exactFuncs_t[TEMP] = pipe_wiggly_T_t;

    1362:       PetscCall(PetscBagGetData(user->bag, (void **)&ctx));
    1363:       id = 2;
    1364:       PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
    1365:       PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    1366:       PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_wiggly_v, 0, NULL));
    1367:       id = 4;
    1368:       PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
    1369:       PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    1370:       PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_wiggly_v, 0, NULL));
    1371:       id = 4;
    1372:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1373:       id = 3;
    1374:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
    1375:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1376:       id = 1;
    1377:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
    1378:       PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
    1379:       break;
    1380:     default:
    1381:       SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
    1382:     }
    1383:     break;
    1384:   default:
    1385:     SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported model type: %s (%d)", modTypes[PetscMin(user->modType, NUM_MOD_TYPES)], user->modType);
    1386:   }
    1387:   /* Setup constants */
    1388:   {
    1389:     Parameter  *param;
    1390:     PetscScalar constants[13];

    1392:     PetscCall(PetscBagGetData(user->bag, (void **)&param));

    1394:     constants[STROUHAL] = param->Strouhal;
    1395:     constants[FROUDE]   = param->Froude;
    1396:     constants[REYNOLDS] = param->Reynolds;
    1397:     constants[PECLET]   = param->Peclet;
    1398:     constants[P_TH]     = param->p_th;
    1399:     constants[MU]       = param->mu;
    1400:     constants[NU]       = param->nu;
    1401:     constants[C_P]      = param->c_p;
    1402:     constants[K]        = param->k;
    1403:     constants[ALPHA]    = param->alpha;
    1404:     constants[T_IN]     = param->T_in;
    1405:     constants[G_DIR]    = param->g_dir;
    1406:     constants[EPSILON]  = param->epsilon;
    1407:     PetscCall(PetscDSSetConstants(ds, 13, constants));
    1408:   }

    1410:   PetscCall(PetscBagGetData(user->bag, (void **)&ctx));
    1411:   PetscCall(PetscDSSetExactSolution(ds, VEL, exactFuncs[VEL], ctx));
    1412:   PetscCall(PetscDSSetExactSolution(ds, PRES, exactFuncs[PRES], ctx));
    1413:   PetscCall(PetscDSSetExactSolution(ds, TEMP, exactFuncs[TEMP], ctx));
    1414:   PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, VEL, exactFuncs_t[VEL], ctx));
    1415:   PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, PRES, exactFuncs_t[PRES], ctx));
    1416:   PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, TEMP, exactFuncs_t[TEMP], ctx));
    1417:   PetscFunctionReturn(PETSC_SUCCESS);
    1418: }

    1420: static PetscErrorCode CreateCellDM(DM dm, AppCtx *user)
    1421: {
    1422:   PetscFE        fe, fediv;
    1423:   DMPolytopeType ct;
    1424:   PetscInt       dim, cStart;
    1425:   PetscBool      simplex;

    1427:   PetscFunctionBeginUser;
    1428:   PetscCall(DMGetDimension(dm, &dim));
    1429:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
    1430:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
    1431:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;

    1433:   PetscCall(DMGetField(dm, VEL, NULL, (PetscObject *)&fe));
    1434:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "div_", PETSC_DEFAULT, &fediv));
    1435:   PetscCall(PetscFECopyQuadrature(fe, fediv));
    1436:   PetscCall(PetscObjectSetName((PetscObject)fediv, "divergence"));

    1438:   PetscCall(DMDestroy(&user->dmCell));
    1439:   PetscCall(DMClone(dm, &user->dmCell));
    1440:   PetscCall(DMSetField(user->dmCell, 0, NULL, (PetscObject)fediv));
    1441:   PetscCall(DMCreateDS(user->dmCell));
    1442:   PetscCall(PetscFEDestroy(&fediv));
    1443:   PetscFunctionReturn(PETSC_SUCCESS);
    1444: }

    1446: static PetscErrorCode GetCellDM(DM dm, AppCtx *user, DM *dmCell)
    1447: {
    1448:   PetscInt cStart, cEnd, cellStart = -1, cellEnd = -1;

    1450:   PetscFunctionBeginUser;
    1451:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
    1452:   if (user->dmCell) PetscCall(DMPlexGetSimplexOrBoxCells(user->dmCell, 0, &cellStart, &cellEnd));
    1453:   if (cStart != cellStart || cEnd != cellEnd) PetscCall(CreateCellDM(dm, user));
    1454:   *dmCell = user->dmCell;
    1455:   PetscFunctionReturn(PETSC_SUCCESS);
    1456: }

    1458: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
    1459: {
    1460:   DM             cdm = dm;
    1461:   PetscFE        fe[3];
    1462:   Parameter     *param;
    1463:   DMPolytopeType ct;
    1464:   PetscInt       dim, cStart;
    1465:   PetscBool      simplex;

    1467:   PetscFunctionBeginUser;
    1468:   PetscCall(DMGetDimension(dm, &dim));
    1469:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
    1470:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
    1471:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
    1472:   /* Create finite element */
    1473:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]));
    1474:   PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));

    1476:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]));
    1477:   PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
    1478:   PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));

    1480:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]));
    1481:   PetscCall(PetscFECopyQuadrature(fe[0], fe[2]));
    1482:   PetscCall(PetscObjectSetName((PetscObject)fe[2], "temperature"));

    1484:   /* Set discretization and boundary conditions for each mesh */
    1485:   PetscCall(DMSetField(dm, VEL, NULL, (PetscObject)fe[VEL]));
    1486:   PetscCall(DMSetField(dm, PRES, NULL, (PetscObject)fe[PRES]));
    1487:   PetscCall(DMSetField(dm, TEMP, NULL, (PetscObject)fe[TEMP]));
    1488:   PetscCall(DMCreateDS(dm));
    1489:   PetscCall(SetupProblem(dm, user));
    1490:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
    1491:   while (cdm) {
    1492:     PetscCall(DMCopyDisc(dm, cdm));
    1493:     PetscCall(DMGetCoarseDM(cdm, &cdm));
    1494:   }
    1495:   PetscCall(PetscFEDestroy(&fe[VEL]));
    1496:   PetscCall(PetscFEDestroy(&fe[PRES]));
    1497:   PetscCall(PetscFEDestroy(&fe[TEMP]));

    1499:   if (user->hasNullSpace) {
    1500:     PetscObject  pressure;
    1501:     MatNullSpace nullspacePres;

    1503:     PetscCall(DMGetField(dm, PRES, NULL, &pressure));
    1504:     PetscCall(MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres));
    1505:     PetscCall(PetscObjectCompose(pressure, "nullspace", (PetscObject)nullspacePres));
    1506:     PetscCall(MatNullSpaceDestroy(&nullspacePres));
    1507:   }
    1508:   PetscFunctionReturn(PETSC_SUCCESS);
    1509: }

    1511: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
    1512: {
    1513:   Vec vec;
    1514:   PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero};

    1516:   PetscFunctionBeginUser;
    1517:   PetscCheck(ofield == PRES, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Nullspace must be for pressure field at index %" PetscInt_FMT ", not %" PetscInt_FMT, PRES, ofield);
    1518:   funcs[nfield] = constant;
    1519:   PetscCall(DMCreateGlobalVector(dm, &vec));
    1520:   PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
    1521:   PetscCall(VecNormalize(vec, NULL));
    1522:   PetscCall(PetscObjectSetName((PetscObject)vec, "Pressure Null Space"));
    1523:   PetscCall(VecViewFromOptions(vec, NULL, "-pressure_nullspace_view"));
    1524:   PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
    1525:   PetscCall(VecDestroy(&vec));
    1526:   PetscFunctionReturn(PETSC_SUCCESS);
    1527: }

    1529: static PetscErrorCode RemoveDiscretePressureNullspace_Private(TS ts, Vec u)
    1530: {
    1531:   DM           dm;
    1532:   AppCtx      *user;
    1533:   MatNullSpace nullsp;

    1535:   PetscFunctionBeginUser;
    1536:   PetscCall(TSGetDM(ts, &dm));
    1537:   PetscCall(DMGetApplicationContext(dm, &user));
    1538:   if (!user->hasNullSpace) PetscFunctionReturn(PETSC_SUCCESS);
    1539:   PetscCall(CreatePressureNullSpace(dm, 1, 1, &nullsp));
    1540:   PetscCall(MatNullSpaceRemove(nullsp, u));
    1541:   PetscCall(MatNullSpaceDestroy(&nullsp));
    1542:   PetscFunctionReturn(PETSC_SUCCESS);
    1543: }

    1545: /* Make the discrete pressure discretely divergence free */
    1546: static PetscErrorCode RemoveDiscretePressureNullspace(TS ts)
    1547: {
    1548:   Vec u;

    1550:   PetscFunctionBeginUser;
    1551:   PetscCall(TSGetSolution(ts, &u));
    1552:   PetscCall(RemoveDiscretePressureNullspace_Private(ts, u));
    1553:   PetscFunctionReturn(PETSC_SUCCESS);
    1554: }

    1556: static void divergence(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 divu[])
    1557: {
    1558:   PetscInt d;

    1560:   divu[0] = 0.;
    1561:   for (d = 0; d < dim; ++d) divu[0] += u_x[d * dim + d];
    1562: }

    1564: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
    1565: {
    1566:   AppCtx   *user;
    1567:   DM        dm;
    1568:   PetscReal t;

    1570:   PetscFunctionBeginUser;
    1571:   PetscCall(TSGetDM(ts, &dm));
    1572:   PetscCall(TSGetTime(ts, &t));
    1573:   PetscCall(DMComputeExactSolution(dm, t, u, NULL));
    1574:   PetscCall(DMGetApplicationContext(dm, &user));
    1575:   PetscCall(RemoveDiscretePressureNullspace_Private(ts, u));
    1576:   PetscFunctionReturn(PETSC_SUCCESS);
    1577: }

    1579: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
    1580: {
    1581:   PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
    1582:   void         *ctxs[3];
    1583:   PetscPointFn *diagnostics[1] = {divergence};
    1584:   DM            dm, dmCell = NULL;
    1585:   PetscDS       ds;
    1586:   Vec           v, divu;
    1587:   PetscReal     ferrors[3], massFlux;
    1588:   PetscInt      f;

    1590:   PetscFunctionBeginUser;
    1591:   PetscCall(TSGetDM(ts, &dm));
    1592:   PetscCall(DMGetDS(dm, &ds));

    1594:   for (f = 0; f < 3; ++f) PetscCall(PetscDSGetExactSolution(ds, f, &exactFuncs[f], &ctxs[f]));
    1595:   PetscCall(DMComputeL2FieldDiff(dm, crtime, exactFuncs, ctxs, u, ferrors));
    1596:   PetscCall(GetCellDM(dm, (AppCtx *)ctx, &dmCell));
    1597:   PetscCall(DMGetGlobalVector(dmCell, &divu));
    1598:   PetscCall(DMProjectField(dmCell, crtime, u, diagnostics, INSERT_VALUES, divu));
    1599:   PetscCall(VecViewFromOptions(divu, NULL, "-divu_vec_view"));
    1600:   PetscCall(VecNorm(divu, NORM_2, &massFlux));
    1601:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g, %2.3g] ||div u||: %2.3g\n", (int)step, (double)crtime, (double)ferrors[0], (double)ferrors[1], (double)ferrors[2], (double)massFlux));

    1603:   PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));

    1605:   PetscCall(DMGetGlobalVector(dm, &v));
    1606:   PetscCall(DMProjectFunction(dm, crtime, exactFuncs, ctxs, INSERT_ALL_VALUES, v));
    1607:   PetscCall(PetscObjectSetName((PetscObject)v, "Exact Solution"));
    1608:   PetscCall(VecViewFromOptions(v, NULL, "-exact_vec_view"));
    1609:   PetscCall(DMRestoreGlobalVector(dm, &v));

    1611:   PetscCall(VecViewFromOptions(divu, NULL, "-div_vec_view"));
    1612:   PetscCall(DMRestoreGlobalVector(dmCell, &divu));
    1613:   PetscFunctionReturn(PETSC_SUCCESS);
    1614: }

    1616: int main(int argc, char **argv)
    1617: {
    1618:   DM        dm;   /* problem definition */
    1619:   TS        ts;   /* timestepper */
    1620:   Vec       u;    /* solution */
    1621:   AppCtx    user; /* user-defined work context */
    1622:   PetscReal t;

    1624:   PetscFunctionBeginUser;
    1625:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    1626:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
    1627:   PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag));
    1628:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
    1629:   PetscCall(SetupParameters(dm, &user));
    1630:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    1631:   PetscCall(TSSetDM(ts, dm));
    1632:   PetscCall(DMSetApplicationContext(dm, &user));
    1633:   /* Setup problem */
    1634:   PetscCall(SetupDiscretization(dm, &user));
    1635:   PetscCall(DMPlexCreateClosureIndex(dm, NULL));

    1637:   PetscCall(DMCreateGlobalVector(dm, &u));
    1638:   PetscCall(PetscObjectSetName((PetscObject)u, "Numerical Solution"));
    1639:   if (user.hasNullSpace) PetscCall(DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace));

    1641:   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user));
    1642:   PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user));
    1643:   PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user));
    1644:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    1645:   PetscCall(TSSetPreStep(ts, RemoveDiscretePressureNullspace));
    1646:   PetscCall(TSSetFromOptions(ts));

    1648:   PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions)); /* Must come after SetFromOptions() */
    1649:   PetscCall(SetInitialConditions(ts, u));
    1650:   PetscCall(TSGetTime(ts, &t));
    1651:   PetscCall(DMSetOutputSequenceNumber(dm, 0, t));
    1652:   PetscCall(DMTSCheckFromOptions(ts, u));
    1653:   PetscCall(TSMonitorSet(ts, MonitorError, &user, NULL));

    1655:   PetscCall(TSSolve(ts, u));
    1656:   PetscCall(DMTSCheckFromOptions(ts, u));

    1658:   PetscCall(VecDestroy(&u));
    1659:   PetscCall(DMDestroy(&user.dmCell));
    1660:   PetscCall(DMDestroy(&dm));
    1661:   PetscCall(TSDestroy(&ts));
    1662:   PetscCall(PetscBagDestroy(&user.bag));
    1663:   PetscCall(PetscFinalize());
    1664:   return 0;
    1665: }

    1667: /*TEST

    1669:   testset:
    1670:     requires: triangle !single
    1671:     args: -dm_plex_separate_marker \
    1672:           -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \
    1673:           -snes_error_if_not_converged -snes_convergence_test correct_pressure \
    1674:           -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
    1675:           -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 \
    1676:           -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
    1677:             -fieldsplit_0_pc_type lu \
    1678:             -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi

    1680:     test:
    1681:       suffix: 2d_tri_p2_p1_p1
    1682:       args: -sol_type quadratic \
    1683:             -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
    1684:             -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1

    1686:     test:
    1687:       # Using -dm_refine 5 -convest_num_refine 2 gives L_2 convergence rate: [0.89, 0.011, 1.0]
    1688:       suffix: 2d_tri_p2_p1_p1_tconv
    1689:       args: -sol_type cubic_trig \
    1690:             -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
    1691:             -ts_max_steps 4 -ts_dt 0.1 -ts_convergence_estimate -convest_num_refine 1

    1693:     test:
    1694:       # Using -dm_refine 3 -convest_num_refine 3 gives L_2 convergence rate: [3.0, 2.5, 1.9]
    1695:       suffix: 2d_tri_p2_p1_p1_sconv
    1696:       args: -sol_type cubic \
    1697:             -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
    1698:             -ts_max_steps 1 -ts_dt 1e-4 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1

    1700:     test:
    1701:       suffix: 2d_tri_p3_p2_p2
    1702:       args: -sol_type cubic \
    1703:             -vel_petscspace_degree 3 -pres_petscspace_degree 2 -temp_petscspace_degree 2 \
    1704:             -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1

    1706:     test:
    1707:       # Using -dm_refine 3 -convest_num_refine 3 gives L_2 convergence rate: [3.0, 2.1, 3.1]
    1708:       suffix: 2d_tri_p2_p1_p1_tg_sconv
    1709:       args: -sol_type taylor_green \
    1710:             -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
    1711:             -ts_max_steps 1 -ts_dt 1e-8 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1

    1713:     test:
    1714:       # Using -dm_refine 3 -convest_num_refine 2 gives L_2 convergence rate: [1.2, 1.5, 1.2]
    1715:       suffix: 2d_tri_p2_p1_p1_tg_tconv
    1716:       args: -sol_type taylor_green \
    1717:             -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
    1718:             -ts_max_steps 4 -ts_dt 0.1 -ts_convergence_estimate -convest_num_refine 1

    1720:   testset:
    1721:     requires: triangle !single
    1722:     args: -dm_plex_separate_marker -mod_type conducting \
    1723:           -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \
    1724:           -snes_error_if_not_converged -snes_max_linear_solve_fail 5 \
    1725:           -ksp_type fgmres -ksp_max_it 2 -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 \
    1726:           -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 \
    1727:           -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
    1728:             -fieldsplit_0_pc_type lu \
    1729:             -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi

    1731:     test:
    1732:       # At this resolution, the rhs is inconsistent on some Newton steps
    1733:       suffix: 2d_tri_p2_p1_p1_conduct
    1734:       args: -sol_type quadratic \
    1735:             -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
    1736:             -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1 \
    1737:             -pc_fieldsplit_schur_precondition full \
    1738:               -fieldsplit_pressure_ksp_max_it 2 -fieldsplit_pressure_pc_type svd

    1740:     test:
    1741:       suffix: 2d_tri_p2_p1_p2_conduct_pipe
    1742:       args: -sol_type pipe \
    1743:             -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 2 \
    1744:             -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1

    1746:     test:
    1747:       suffix: 2d_tri_p2_p1_p2_conduct_pipe_wiggly_sconv
    1748:       args: -sol_type pipe_wiggly -Fr 1e10 \
    1749:             -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 2 \
    1750:             -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
    1751:             -ts_max_steps 1 -ts_dt 1e10 \
    1752:             -ksp_atol 1e-12 -ksp_max_it 300 \
    1753:               -fieldsplit_pressure_ksp_atol 1e-14

    1755: TEST*/