Actual source code: ex70.c

      1: static char help[] = "------------------------------------------------------------------------------------------------------------------------------ \n\
      2:   Solves the time-dependent incompressible, variable viscosity Stokes equation in 2D driven by buoyancy variations. \n\
      3:   Time-dependence is introduced by evolving the density (rho) and viscosity (eta) according to \n\
      4:     D \\rho / Dt = 0    and    D \\eta / Dt = 0 \n\
      5:   The Stokes problem is discretized using Q1-Q1 finite elements, stabilized with Bochev's polynomial projection method. \n\
      6:   The hyperbolic evolution equation for density is discretized using a variant of the Particle-In-Cell (PIC) method. \n\
      7:   The DMDA object is used to define the FE problem, whilst DMSwarm provides support for the PIC method. \n\
      8:   Material points (particles) store density and viscosity. The particles are advected with the fluid velocity using RK1. \n\
      9:   At each time step, the value of density and viscosity stored on each particle are projected into a Q1 function space \n\
     10:   and then interpolated onto the Gauss quadrature points. \n\
     11:   The model problem defined in this example is the iso-viscous Rayleigh-Taylor instability (case 1a) from: \n\
     12:     \"A comparison of methods for the modeling of thermochemical convection\" \n\
     13:     P.E. van Keken, S.D. King, H. Schmeling, U.R. Christensen, D. Neumeister and M.-P. Doin, \n\
     14:     Journal of Geophysical Research, vol 102 (B10), 477--499 (1997) \n\
     15:   Note that whilst the model problem defined is for an iso-viscous, the implementation in this example supports \n\
     16:   variable viscosity formulations. \n\
     17:   This example is based on src/ksp/ksp/tutorials/ex43.c \n\
     18:   Options: \n\
     19:     -mx        : Number of elements in the x-direction \n\
     20:     -my        : Number of elements in the y-direction \n\
     21:     -mxy       : Number of elements in the x- and y-directions \n\
     22:     -nt        : Number of time steps \n\
     23:     -dump_freq : Frequency of output file creation \n\
     24:     -ppcell    : Number of times the reference cell is sub-divided \n\
     25:     -randomize_coords : Apply a random shift to each particle coordinate in the range [-fac*dh,0.fac*dh] \n\
     26:     -randomize_fac    : Set the scaling factor for the random shift (default = 0.25)\n";

     28: /* Contributed by Dave May */

     30: #include <petscksp.h>
     31: #include <petscdm.h>
     32: #include <petscdmda.h>
     33: #include <petscdmswarm.h>
     34: #include <petsc/private/dmimpl.h>

     36: static PetscErrorCode DMDAApplyBoundaryConditions(DM, Mat, Vec);

     38: #define NSD          2 /* number of spatial dimensions */
     39: #define NODES_PER_EL 4 /* nodes per element */
     40: #define U_DOFS       2 /* degrees of freedom per velocity node */
     41: #define P_DOFS       1 /* degrees of freedom per pressure node */
     42: #define GAUSS_POINTS 4

     44: static void EvaluateBasis_Q1(PetscScalar _xi[], PetscScalar N[])
     45: {
     46:   PetscScalar xi  = _xi[0];
     47:   PetscScalar eta = _xi[1];

     49:   N[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
     50:   N[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
     51:   N[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
     52:   N[3] = 0.25 * (1.0 - xi) * (1.0 + eta);
     53: }

     55: static void EvaluateBasisDerivatives_Q1(PetscScalar _xi[], PetscScalar dN[][NODES_PER_EL])
     56: {
     57:   PetscScalar xi  = _xi[0];
     58:   PetscScalar eta = _xi[1];

     60:   dN[0][0] = -0.25 * (1.0 - eta);
     61:   dN[0][1] = 0.25 * (1.0 - eta);
     62:   dN[0][2] = 0.25 * (1.0 + eta);
     63:   dN[0][3] = -0.25 * (1.0 + eta);

     65:   dN[1][0] = -0.25 * (1.0 - xi);
     66:   dN[1][1] = -0.25 * (1.0 + xi);
     67:   dN[1][2] = 0.25 * (1.0 + xi);
     68:   dN[1][3] = 0.25 * (1.0 - xi);
     69: }

     71: static void EvaluateDerivatives(PetscScalar dN[][NODES_PER_EL], PetscScalar dNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
     72: {
     73:   PetscScalar J00, J01, J10, J11, J;
     74:   PetscScalar iJ00, iJ01, iJ10, iJ11;
     75:   PetscInt    i;

     77:   J00 = J01 = J10 = J11 = 0.0;
     78:   for (i = 0; i < NODES_PER_EL; i++) {
     79:     PetscScalar cx = coords[2 * i];
     80:     PetscScalar cy = coords[2 * i + 1];

     82:     J00 += dN[0][i] * cx; /* J_xx = dx/dxi */
     83:     J01 += dN[0][i] * cy; /* J_xy = dy/dxi */
     84:     J10 += dN[1][i] * cx; /* J_yx = dx/deta */
     85:     J11 += dN[1][i] * cy; /* J_yy = dy/deta */
     86:   }
     87:   J = (J00 * J11) - (J01 * J10);

     89:   iJ00 = J11 / J;
     90:   iJ01 = -J01 / J;
     91:   iJ10 = -J10 / J;
     92:   iJ11 = J00 / J;

     94:   for (i = 0; i < NODES_PER_EL; i++) {
     95:     dNx[0][i] = dN[0][i] * iJ00 + dN[1][i] * iJ01;
     96:     dNx[1][i] = dN[0][i] * iJ10 + dN[1][i] * iJ11;
     97:   }

     99:   if (det_J) *det_J = J;
    100: }

    102: static void CreateGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
    103: {
    104:   *ngp         = 4;
    105:   gp_xi[0][0]  = -0.57735026919;
    106:   gp_xi[0][1]  = -0.57735026919;
    107:   gp_xi[1][0]  = -0.57735026919;
    108:   gp_xi[1][1]  = 0.57735026919;
    109:   gp_xi[2][0]  = 0.57735026919;
    110:   gp_xi[2][1]  = 0.57735026919;
    111:   gp_xi[3][0]  = 0.57735026919;
    112:   gp_xi[3][1]  = -0.57735026919;
    113:   gp_weight[0] = 1.0;
    114:   gp_weight[1] = 1.0;
    115:   gp_weight[2] = 1.0;
    116:   gp_weight[3] = 1.0;
    117: }

    119: static PetscErrorCode DMDAGetElementEqnums_up(const PetscInt element[], PetscInt s_u[], PetscInt s_p[])
    120: {
    121:   PetscFunctionBeginUser;
    122:   for (PetscInt i = 0; i < NODES_PER_EL; i++) {
    123:     /* velocity */
    124:     s_u[NSD * i + 0] = 3 * element[i];
    125:     s_u[NSD * i + 1] = 3 * element[i] + 1;
    126:     /* pressure */
    127:     s_p[i] = 3 * element[i] + 2;
    128:   }
    129:   PetscFunctionReturn(PETSC_SUCCESS);
    130: }

    132: static PetscInt map_wIwDI_uJuDJ(PetscInt wi, PetscInt wd, PetscInt w_NPE, PetscInt w_dof, PetscInt ui, PetscInt ud, PetscInt u_NPE, PetscInt u_dof)
    133: {
    134:   PetscInt ij, r, c, nc;

    136:   nc = u_NPE * u_dof;
    137:   r  = w_dof * wi + wd;
    138:   c  = u_dof * ui + ud;
    139:   ij = r * nc + c;
    140:   return ij;
    141: }

    143: static void BForm_DivT(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
    144: {
    145:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
    146:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    147:   PetscScalar J_p, tildeD[3];
    148:   PetscScalar B[3][U_DOFS * NODES_PER_EL];
    149:   PetscInt    p, i, j, k, ngp;

    151:   /* define quadrature rule */
    152:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

    154:   /* evaluate bilinear form */
    155:   for (p = 0; p < ngp; p++) {
    156:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
    157:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);

    159:     for (i = 0; i < NODES_PER_EL; i++) {
    160:       PetscScalar d_dx_i = GNx_p[0][i];
    161:       PetscScalar d_dy_i = GNx_p[1][i];

    163:       B[0][2 * i]     = d_dx_i;
    164:       B[0][2 * i + 1] = 0.0;
    165:       B[1][2 * i]     = 0.0;
    166:       B[1][2 * i + 1] = d_dy_i;
    167:       B[2][2 * i]     = d_dy_i;
    168:       B[2][2 * i + 1] = d_dx_i;
    169:     }

    171:     tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
    172:     tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
    173:     tildeD[2] = gp_weight[p] * J_p * eta[p];

    175:     /* form Bt tildeD B */
    176:     /*
    177:     Ke_ij = Bt_ik . D_kl . B_lj
    178:           = B_ki . D_kl . B_lj
    179:           = B_ki . D_kk . B_kj
    180:     */
    181:     for (i = 0; i < 8; i++) {
    182:       for (j = 0; j < 8; j++) {
    183:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
    184:           Ke[i + 8 * j] += B[k][i] * tildeD[k] * B[k][j];
    185:         }
    186:       }
    187:     }
    188:   }
    189: }

    191: static void BForm_Grad(PetscScalar Ke[], PetscScalar coords[])
    192: {
    193:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
    194:   PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    195:   PetscScalar J_p, fac;
    196:   PetscInt    p, i, j, di, ngp;

    198:   /* define quadrature rule */
    199:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

    201:   /* evaluate bilinear form */
    202:   for (p = 0; p < ngp; p++) {
    203:     EvaluateBasis_Q1(gp_xi[p], Ni_p);
    204:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
    205:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
    206:     fac = gp_weight[p] * J_p;

    208:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
    209:       for (di = 0; di < NSD; di++) {     /* u dofs */
    210:         for (j = 0; j < 4; j++) {        /* p nodes, p dofs = 1 (ie no loop) */
    211:           PetscInt IJ;
    212:           IJ = map_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);

    214:           Ke[IJ] -= GNx_p[di][i] * Ni_p[j] * fac;
    215:         }
    216:       }
    217:     }
    218:   }
    219: }

    221: static void BForm_Div(PetscScalar De[], PetscScalar coords[])
    222: {
    223:   PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
    224:   PetscInt    i, j, nr_g, nc_g;

    226:   PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
    227:   BForm_Grad(Ge, coords);

    229:   nr_g = U_DOFS * NODES_PER_EL;
    230:   nc_g = P_DOFS * NODES_PER_EL;

    232:   for (i = 0; i < nr_g; i++) {
    233:     for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
    234:   }
    235: }

    237: static void BForm_Stabilisation(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
    238: {
    239:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
    240:   PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    241:   PetscScalar J_p, fac, eta_avg;
    242:   PetscInt    p, i, j, ngp;

    244:   /* define quadrature rule */
    245:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

    247:   /* evaluate bilinear form */
    248:   for (p = 0; p < ngp; p++) {
    249:     EvaluateBasis_Q1(gp_xi[p], Ni_p);
    250:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
    251:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
    252:     fac = gp_weight[p] * J_p;

    254:     for (i = 0; i < NODES_PER_EL; i++) {
    255:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * (Ni_p[i] * Ni_p[j] - 0.0625);
    256:     }
    257:   }

    259:   /* scale */
    260:   eta_avg = 0.0;
    261:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    262:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
    263:   fac     = 1.0 / eta_avg;
    264:   for (i = 0; i < NODES_PER_EL; i++) {
    265:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
    266:   }
    267: }

    269: static void BForm_ScaledMassMatrix(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
    270: {
    271:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
    272:   PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    273:   PetscScalar J_p, fac, eta_avg;
    274:   PetscInt    p, i, j, ngp;

    276:   /* define quadrature rule */
    277:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

    279:   /* evaluate bilinear form */
    280:   for (p = 0; p < ngp; p++) {
    281:     EvaluateBasis_Q1(gp_xi[p], Ni_p);
    282:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
    283:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
    284:     fac = gp_weight[p] * J_p;

    286:     for (i = 0; i < NODES_PER_EL; i++) {
    287:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * Ni_p[i] * Ni_p[j];
    288:     }
    289:   }

    291:   /* scale */
    292:   eta_avg = 0.0;
    293:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    294:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
    295:   fac     = 1.0 / eta_avg;
    296:   for (i = 0; i < NODES_PER_EL; i++) {
    297:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] *= fac;
    298:   }
    299: }

    301: static void LForm_MomentumRHS(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
    302: {
    303:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
    304:   PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    305:   PetscScalar J_p, fac;
    306:   PetscInt    p, i, ngp;

    308:   /* define quadrature rule */
    309:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

    311:   /* evaluate linear form */
    312:   for (p = 0; p < ngp; p++) {
    313:     EvaluateBasis_Q1(gp_xi[p], Ni_p);
    314:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
    315:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
    316:     fac = gp_weight[p] * J_p;

    318:     for (i = 0; i < NODES_PER_EL; i++) {
    319:       Fe[NSD * i] = 0.0;
    320:       Fe[NSD * i + 1] -= fac * Ni_p[i] * fy[p];
    321:     }
    322:   }
    323: }

    325: static PetscErrorCode GetElementCoords(const PetscScalar _coords[], const PetscInt e2n[], PetscScalar el_coords[])
    326: {
    327:   PetscFunctionBeginUser;
    328:   /* get coords for the element */
    329:   for (PetscInt i = 0; i < 4; i++) {
    330:     for (PetscInt d = 0; d < NSD; d++) el_coords[NSD * i + d] = _coords[NSD * e2n[i] + d];
    331:   }
    332:   PetscFunctionReturn(PETSC_SUCCESS);
    333: }

    335: static PetscErrorCode AssembleStokes_A(Mat A, DM stokes_da, DM quadrature)
    336: {
    337:   DM                 cda;
    338:   Vec                coords;
    339:   const PetscScalar *_coords;
    340:   PetscInt           u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
    341:   PetscInt           p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
    342:   PetscInt           nel, npe, eidx;
    343:   const PetscInt    *element_list;
    344:   PetscScalar        Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
    345:   PetscScalar        Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
    346:   PetscScalar        De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
    347:   PetscScalar        Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
    348:   PetscScalar        el_coords[NODES_PER_EL * NSD];
    349:   PetscScalar       *q_eta, *prop_eta;

    351:   PetscFunctionBeginUser;
    352:   PetscCall(MatZeroEntries(A));
    353:   /* setup for coords */
    354:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    355:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    356:   PetscCall(VecGetArrayRead(coords, &_coords));

    358:   /* setup for coefficients */
    359:   PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));

    361:   PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
    362:   for (eidx = 0; eidx < nel; eidx++) {
    363:     const PetscInt *element = &element_list[npe * eidx];

    365:     /* get coords for the element */
    366:     PetscCall(GetElementCoords(_coords, element, el_coords));

    368:     /* get coefficients for the element */
    369:     prop_eta = &q_eta[GAUSS_POINTS * eidx];

    371:     /* initialise element stiffness matrix */
    372:     PetscCall(PetscMemzero(Ae, sizeof(Ae)));
    373:     PetscCall(PetscMemzero(Ge, sizeof(Ge)));
    374:     PetscCall(PetscMemzero(De, sizeof(De)));
    375:     PetscCall(PetscMemzero(Ce, sizeof(Ce)));

    377:     /* form element stiffness matrix */
    378:     BForm_DivT(Ae, el_coords, prop_eta);
    379:     BForm_Grad(Ge, el_coords);
    380:     BForm_Div(De, el_coords);
    381:     BForm_Stabilisation(Ce, el_coords, prop_eta);

    383:     /* insert element matrix into global matrix */
    384:     PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
    385:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
    386:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
    387:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
    388:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
    389:   }
    390:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    391:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    393:   PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
    394:   PetscCall(VecRestoreArrayRead(coords, &_coords));
    395:   PetscFunctionReturn(PETSC_SUCCESS);
    396: }

    398: static PetscErrorCode AssembleStokes_PC(Mat A, DM stokes_da, DM quadrature)
    399: {
    400:   DM                 cda;
    401:   Vec                coords;
    402:   const PetscScalar *_coords;
    403:   PetscInt           u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
    404:   PetscInt           p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
    405:   PetscInt           nel, npe, eidx;
    406:   const PetscInt    *element_list;
    407:   PetscScalar        Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
    408:   PetscScalar        Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
    409:   PetscScalar        De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
    410:   PetscScalar        Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
    411:   PetscScalar        el_coords[NODES_PER_EL * NSD];
    412:   PetscScalar       *q_eta, *prop_eta;

    414:   PetscFunctionBeginUser;
    415:   PetscCall(MatZeroEntries(A));
    416:   /* setup for coords */
    417:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    418:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    419:   PetscCall(VecGetArrayRead(coords, &_coords));

    421:   /* setup for coefficients */
    422:   PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));

    424:   PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
    425:   for (eidx = 0; eidx < nel; eidx++) {
    426:     const PetscInt *element = &element_list[npe * eidx];

    428:     /* get coords for the element */
    429:     PetscCall(GetElementCoords(_coords, element, el_coords));

    431:     /* get coefficients for the element */
    432:     prop_eta = &q_eta[GAUSS_POINTS * eidx];

    434:     /* initialise element stiffness matrix */
    435:     PetscCall(PetscMemzero(Ae, sizeof(Ae)));
    436:     PetscCall(PetscMemzero(Ge, sizeof(Ge)));
    437:     PetscCall(PetscMemzero(De, sizeof(De)));
    438:     PetscCall(PetscMemzero(Ce, sizeof(Ce)));

    440:     /* form element stiffness matrix */
    441:     BForm_DivT(Ae, el_coords, prop_eta);
    442:     BForm_Grad(Ge, el_coords);
    443:     BForm_ScaledMassMatrix(Ce, el_coords, prop_eta);

    445:     /* insert element matrix into global matrix */
    446:     PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
    447:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
    448:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
    449:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
    450:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
    451:   }
    452:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    453:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    455:   PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
    456:   PetscCall(VecRestoreArrayRead(coords, &_coords));
    457:   PetscFunctionReturn(PETSC_SUCCESS);
    458: }

    460: static PetscErrorCode AssembleStokes_RHS(Vec F, DM stokes_da, DM quadrature)
    461: {
    462:   DM                 cda;
    463:   Vec                coords;
    464:   const PetscScalar *_coords;
    465:   PetscInt           u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
    466:   PetscInt           p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
    467:   PetscInt           nel, npe, eidx, i;
    468:   const PetscInt    *element_list;
    469:   PetscScalar        Fe[NODES_PER_EL * U_DOFS];
    470:   PetscScalar        He[NODES_PER_EL * P_DOFS];
    471:   PetscScalar        el_coords[NODES_PER_EL * NSD];
    472:   PetscScalar       *q_rhs, *prop_fy;
    473:   Vec                local_F;
    474:   PetscScalar       *LA_F;

    476:   PetscFunctionBeginUser;
    477:   PetscCall(VecZeroEntries(F));
    478:   /* setup for coords */
    479:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    480:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    481:   PetscCall(VecGetArrayRead(coords, &_coords));

    483:   /* setup for coefficients */
    484:   PetscCall(DMSwarmGetField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));

    486:   /* get access to the vector */
    487:   PetscCall(DMGetLocalVector(stokes_da, &local_F));
    488:   PetscCall(VecZeroEntries(local_F));
    489:   PetscCall(VecGetArray(local_F, &LA_F));

    491:   PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
    492:   for (eidx = 0; eidx < nel; eidx++) {
    493:     const PetscInt *element = &element_list[npe * eidx];

    495:     /* get coords for the element */
    496:     PetscCall(GetElementCoords(_coords, element, el_coords));

    498:     /* get coefficients for the element */
    499:     prop_fy = &q_rhs[GAUSS_POINTS * eidx];

    501:     /* initialise element stiffness matrix */
    502:     PetscCall(PetscMemzero(Fe, sizeof(Fe)));
    503:     PetscCall(PetscMemzero(He, sizeof(He)));

    505:     /* form element stiffness matrix */
    506:     LForm_MomentumRHS(Fe, el_coords, NULL, prop_fy);

    508:     /* insert element matrix into global matrix */
    509:     PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));

    511:     for (i = 0; i < NODES_PER_EL * U_DOFS; i++) LA_F[u_eqn[i]] += Fe[i];
    512:   }
    513:   PetscCall(DMSwarmRestoreField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
    514:   PetscCall(VecRestoreArrayRead(coords, &_coords));

    516:   PetscCall(VecRestoreArray(local_F, &LA_F));
    517:   PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
    518:   PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
    519:   PetscCall(DMRestoreLocalVector(stokes_da, &local_F));
    520:   PetscFunctionReturn(PETSC_SUCCESS);
    521: }

    523: PetscErrorCode DMSwarmPICInsertPointsCellwise(DM dm, DM dmc, PetscInt e, PetscInt npoints, PetscReal xi[], PetscBool proximity_initialization)
    524: {
    525:   DMSwarmCellDM      celldm;
    526:   PetscInt           dim, nel, npe, q, k, d, ncurr, Nfc;
    527:   const PetscInt    *element_list;
    528:   Vec                coor;
    529:   const PetscScalar *_coor;
    530:   PetscReal        **basis, *elcoor, *xp;
    531:   PetscReal         *swarm_coor;
    532:   PetscInt          *swarm_cellid;
    533:   const char       **coordFields, *cellid;

    535:   PetscFunctionBeginUser;
    536:   PetscCall(DMGetDimension(dm, &dim));
    537:   PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list));

    539:   PetscCall(PetscMalloc1(dim * npoints, &xp));
    540:   PetscCall(PetscMalloc1(dim * npe, &elcoor));
    541:   PetscCall(PetscMalloc1(npoints, &basis));
    542:   for (q = 0; q < npoints; q++) {
    543:     PetscCall(PetscMalloc1(npe, &basis[q]));

    545:     switch (dim) {
    546:     case 1:
    547:       basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]);
    548:       basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]);
    549:       break;
    550:     case 2:
    551:       basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
    552:       basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
    553:       basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
    554:       basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
    555:       break;

    557:     case 3:
    558:       basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
    559:       basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
    560:       basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
    561:       basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
    562:       basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
    563:       basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
    564:       basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
    565:       basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
    566:       break;
    567:     }
    568:   }

    570:   PetscCall(DMGetCoordinatesLocal(dmc, &coor));
    571:   PetscCall(VecGetArrayRead(coor, &_coor));
    572:   /* compute and store the coordinates for the new points */
    573:   {
    574:     const PetscInt *element = &element_list[npe * e];

    576:     for (k = 0; k < npe; k++) {
    577:       for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
    578:     }
    579:     for (q = 0; q < npoints; q++) {
    580:       for (d = 0; d < dim; d++) xp[dim * q + d] = 0.0;
    581:       for (k = 0; k < npe; k++) {
    582:         for (d = 0; d < dim; d++) xp[dim * q + d] += basis[q][k] * elcoor[dim * k + d];
    583:       }
    584:     }
    585:   }
    586:   PetscCall(VecRestoreArrayRead(coor, &_coor));
    587:   PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list));

    589:   PetscCall(DMSwarmGetLocalSize(dm, &ncurr));
    590:   PetscCall(DMSwarmAddNPoints(dm, npoints));
    591:   PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
    592:   PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
    593:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
    594:   PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);

    596:   if (proximity_initialization) {
    597:     PetscInt  *nnlist;
    598:     PetscReal *coor_q, *coor_qn;
    599:     PetscInt   npoints_e, *plist_e;

    601:     PetscCall(DMSwarmSortGetPointsPerCell(dm, e, &npoints_e, &plist_e));

    603:     PetscCall(PetscMalloc1(npoints, &nnlist));
    604:     /* find nearest neighbour points in this cell */
    605:     PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
    606:     PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
    607:     for (q = 0; q < npoints; q++) {
    608:       PetscInt  qn, nearest_neighbour = -1;
    609:       PetscReal sep, min_sep          = PETSC_MAX_REAL;

    611:       coor_q = &xp[dim * q];
    612:       for (qn = 0; qn < npoints_e; qn++) {
    613:         coor_qn = &swarm_coor[dim * plist_e[qn]];
    614:         sep     = 0.0;
    615:         for (d = 0; d < dim; d++) sep += (coor_q[d] - coor_qn[d]) * (coor_q[d] - coor_qn[d]);
    616:         if (sep < min_sep) {
    617:           nearest_neighbour = plist_e[qn];
    618:           min_sep           = sep;
    619:         }
    620:       }
    621:       PetscCheck(nearest_neighbour != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell %" PetscInt_FMT " is empty - cannot initialize using nearest neighbours", e);
    622:       nnlist[q] = nearest_neighbour;
    623:     }
    624:     PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
    625:     PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));

    627:     /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
    628:     for (q = 0; q < npoints; q++) PetscCall(DMSwarmCopyPoint(dm, nnlist[q], ncurr + q));
    629:     PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
    630:     PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
    631:     for (q = 0; q < npoints; q++) {
    632:       /* set the coordinates */
    633:       for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
    634:       /* set the cell index */
    635:       swarm_cellid[ncurr + q] = e;
    636:     }
    637:     PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
    638:     PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));

    640:     PetscCall(DMSwarmSortRestorePointsPerCell(dm, e, &npoints_e, &plist_e));
    641:     PetscCall(PetscFree(nnlist));
    642:   } else {
    643:     PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
    644:     PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
    645:     for (q = 0; q < npoints; q++) {
    646:       /* set the coordinates */
    647:       for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
    648:       /* set the cell index */
    649:       swarm_cellid[ncurr + q] = e;
    650:     }
    651:     PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
    652:     PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
    653:   }

    655:   PetscCall(PetscFree(xp));
    656:   PetscCall(PetscFree(elcoor));
    657:   for (q = 0; q < npoints; q++) PetscCall(PetscFree(basis[q]));
    658:   PetscCall(PetscFree(basis));
    659:   PetscFunctionReturn(PETSC_SUCCESS);
    660: }

    662: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp, DM dm_mpoint)
    663: {
    664:   PetscInt         _npe, _nel, e, nel;
    665:   const PetscInt  *element;
    666:   DM               dmc;
    667:   PetscQuadrature  quadrature;
    668:   const PetscReal *xi;
    669:   PetscInt         npoints_q, cnt, cnt_g;

    671:   PetscFunctionBeginUser;
    672:   PetscCall(DMDAGetElements(dm_vp, &_nel, &_npe, &element));
    673:   nel = _nel;
    674:   PetscCall(DMDARestoreElements(dm_vp, &_nel, &_npe, &element));

    676:   PetscCall(PetscDTGaussTensorQuadrature(2, 1, 4, -1.0, 1.0, &quadrature));
    677:   PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xi, NULL));
    678:   PetscCall(DMSwarmGetCellDM(dm_mpoint, &dmc));

    680:   PetscCall(DMSwarmSortGetAccess(dm_mpoint));

    682:   cnt = 0;
    683:   for (e = 0; e < nel; e++) {
    684:     PetscInt npoints_per_cell;

    686:     PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint, e, &npoints_per_cell));

    688:     if (npoints_per_cell < 12) {
    689:       PetscCall(DMSwarmPICInsertPointsCellwise(dm_mpoint, dm_vp, e, npoints_q, (PetscReal *)xi, PETSC_TRUE));
    690:       cnt++;
    691:     }
    692:   }
    693:   PetscCallMPI(MPIU_Allreduce(&cnt, &cnt_g, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
    694:   if (cnt_g > 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... ....pop cont: adjusted %" PetscInt_FMT " cells\n", cnt_g));

    696:   PetscCall(DMSwarmSortRestoreAccess(dm_mpoint));
    697:   PetscCall(PetscQuadratureDestroy(&quadrature));
    698:   PetscFunctionReturn(PETSC_SUCCESS);
    699: }

    701: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp, Vec vp, PetscReal dt, DM dm_mpoint)
    702: {
    703:   DMSwarmCellDM      celldm;
    704:   Vec                vp_l, coor_l;
    705:   const PetscScalar *LA_vp;
    706:   PetscInt           i, p, e, npoints, nel, npe, Nfc;
    707:   PetscInt          *mpfield_cell;
    708:   PetscReal         *mpfield_coor;
    709:   const PetscInt    *element_list;
    710:   const PetscInt    *element;
    711:   PetscScalar        xi_p[NSD], Ni[NODES_PER_EL];
    712:   const PetscScalar *LA_coor;
    713:   PetscScalar        dx[NSD];
    714:   const char       **coordFields, *cellid;

    716:   PetscFunctionBeginUser;
    717:   PetscCall(DMGetCoordinatesLocal(dm_vp, &coor_l));
    718:   PetscCall(VecGetArrayRead(coor_l, &LA_coor));

    720:   PetscCall(DMGetLocalVector(dm_vp, &vp_l));
    721:   PetscCall(DMGlobalToLocalBegin(dm_vp, vp, INSERT_VALUES, vp_l));
    722:   PetscCall(DMGlobalToLocalEnd(dm_vp, vp, INSERT_VALUES, vp_l));
    723:   PetscCall(VecGetArrayRead(vp_l, &LA_vp));

    725:   PetscCall(DMDAGetElements(dm_vp, &nel, &npe, &element_list));
    726:   PetscCall(DMSwarmGetLocalSize(dm_mpoint, &npoints));
    727:   PetscCall(DMSwarmGetCellDMActive(dm_mpoint, &celldm));
    728:   PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
    729:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
    730:   PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm_mpoint), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
    731:   PetscCall(DMSwarmGetField(dm_mpoint, coordFields[0], NULL, NULL, (void **)&mpfield_coor));
    732:   PetscCall(DMSwarmGetField(dm_mpoint, cellid, NULL, NULL, (void **)&mpfield_cell));
    733:   for (p = 0; p < npoints; p++) {
    734:     PetscReal         *coor_p;
    735:     PetscScalar        vel_n[NSD * NODES_PER_EL], vel_p[NSD];
    736:     const PetscScalar *x0;
    737:     const PetscScalar *x2;

    739:     e       = mpfield_cell[p];
    740:     coor_p  = &mpfield_coor[NSD * p];
    741:     element = &element_list[NODES_PER_EL * e];

    743:     /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
    744:     x0 = &LA_coor[NSD * element[0]];
    745:     x2 = &LA_coor[NSD * element[2]];

    747:     dx[0] = x2[0] - x0[0];
    748:     dx[1] = x2[1] - x0[1];

    750:     xi_p[0] = 2.0 * (coor_p[0] - x0[0]) / dx[0] - 1.0;
    751:     xi_p[1] = 2.0 * (coor_p[1] - x0[1]) / dx[1] - 1.0;
    752:     PetscCheck(PetscRealPart(xi_p[0]) >= -1.0 - PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (xi) too small %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[0]), e);
    753:     PetscCheck(PetscRealPart(xi_p[0]) <= 1.0 + PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (xi) too large %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[0]), e);
    754:     PetscCheck(PetscRealPart(xi_p[1]) >= -1.0 - PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (eta) too small %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[1]), e);
    755:     PetscCheck(PetscRealPart(xi_p[1]) <= 1.0 + PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (eta) too large %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[1]), e);

    757:     /* evaluate basis functions */
    758:     EvaluateBasis_Q1(xi_p, Ni);

    760:     /* get cell nodal velocities */
    761:     for (i = 0; i < NODES_PER_EL; i++) {
    762:       PetscInt nid;

    764:       nid                = element[i];
    765:       vel_n[NSD * i + 0] = LA_vp[(NSD + 1) * nid + 0];
    766:       vel_n[NSD * i + 1] = LA_vp[(NSD + 1) * nid + 1];
    767:     }

    769:     /* interpolate velocity */
    770:     vel_p[0] = vel_p[1] = 0.0;
    771:     for (i = 0; i < NODES_PER_EL; i++) {
    772:       vel_p[0] += Ni[i] * vel_n[NSD * i + 0];
    773:       vel_p[1] += Ni[i] * vel_n[NSD * i + 1];
    774:     }

    776:     coor_p[0] += dt * PetscRealPart(vel_p[0]);
    777:     coor_p[1] += dt * PetscRealPart(vel_p[1]);
    778:   }

    780:   PetscCall(DMSwarmRestoreField(dm_mpoint, cellid, NULL, NULL, (void **)&mpfield_cell));
    781:   PetscCall(DMSwarmRestoreField(dm_mpoint, coordFields[0], NULL, NULL, (void **)&mpfield_coor));
    782:   PetscCall(DMDARestoreElements(dm_vp, &nel, &npe, &element_list));
    783:   PetscCall(VecRestoreArrayRead(vp_l, &LA_vp));
    784:   PetscCall(DMRestoreLocalVector(dm_vp, &vp_l));
    785:   PetscCall(VecRestoreArrayRead(coor_l, &LA_coor));
    786:   PetscFunctionReturn(PETSC_SUCCESS);
    787: }

    789: PetscErrorCode MaterialPoint_Interpolate(DM dm, Vec eta_v, Vec rho_v, DM dm_quadrature)
    790: {
    791:   Vec             eta_l, rho_l;
    792:   PetscScalar    *_eta_l, *_rho_l;
    793:   PetscInt        nqp, npe, nel;
    794:   PetscScalar     qp_xi[GAUSS_POINTS][NSD];
    795:   PetscScalar     qp_weight[GAUSS_POINTS];
    796:   PetscInt        q, k, e;
    797:   PetscScalar     Ni[GAUSS_POINTS][NODES_PER_EL];
    798:   const PetscInt *element_list;
    799:   PetscReal      *q_eta, *q_rhs;

    801:   PetscFunctionBeginUser;
    802:   /* define quadrature rule */
    803:   CreateGaussQuadrature(&nqp, qp_xi, qp_weight);
    804:   for (q = 0; q < nqp; q++) EvaluateBasis_Q1(qp_xi[q], Ni[q]);

    806:   PetscCall(DMGetLocalVector(dm, &eta_l));
    807:   PetscCall(DMGetLocalVector(dm, &rho_l));

    809:   PetscCall(DMGlobalToLocalBegin(dm, eta_v, INSERT_VALUES, eta_l));
    810:   PetscCall(DMGlobalToLocalEnd(dm, eta_v, INSERT_VALUES, eta_l));
    811:   PetscCall(DMGlobalToLocalBegin(dm, rho_v, INSERT_VALUES, rho_l));
    812:   PetscCall(DMGlobalToLocalEnd(dm, rho_v, INSERT_VALUES, rho_l));

    814:   PetscCall(VecGetArray(eta_l, &_eta_l));
    815:   PetscCall(VecGetArray(rho_l, &_rho_l));

    817:   PetscCall(DMSwarmGetField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
    818:   PetscCall(DMSwarmGetField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));

    820:   PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list));
    821:   for (e = 0; e < nel; e++) {
    822:     PetscScalar     eta_field_e[NODES_PER_EL];
    823:     PetscScalar     rho_field_e[NODES_PER_EL];
    824:     const PetscInt *element = &element_list[4 * e];

    826:     for (k = 0; k < NODES_PER_EL; k++) {
    827:       eta_field_e[k] = _eta_l[element[k]];
    828:       rho_field_e[k] = _rho_l[element[k]];
    829:     }

    831:     for (q = 0; q < nqp; q++) {
    832:       PetscScalar eta_q, rho_q;

    834:       eta_q = rho_q = 0.0;
    835:       for (k = 0; k < NODES_PER_EL; k++) {
    836:         eta_q += Ni[q][k] * eta_field_e[k];
    837:         rho_q += Ni[q][k] * rho_field_e[k];
    838:       }

    840:       q_eta[nqp * e + q] = PetscRealPart(eta_q);
    841:       q_rhs[nqp * e + q] = PetscRealPart(rho_q);
    842:     }
    843:   }
    844:   PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list));

    846:   PetscCall(DMSwarmRestoreField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
    847:   PetscCall(DMSwarmRestoreField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));

    849:   PetscCall(VecRestoreArray(rho_l, &_rho_l));
    850:   PetscCall(VecRestoreArray(eta_l, &_eta_l));
    851:   PetscCall(DMRestoreLocalVector(dm, &rho_l));
    852:   PetscCall(DMRestoreLocalVector(dm, &eta_l));
    853:   PetscFunctionReturn(PETSC_SUCCESS);
    854: }

    856: static PetscErrorCode SolveTimeDepStokes(PetscInt mx, PetscInt my)
    857: {
    858:   DM              dm_stokes, dm_coeff;
    859:   PetscInt        u_dof, p_dof, dof, stencil_width;
    860:   Mat             A, B;
    861:   PetscInt        nel_local;
    862:   Vec             eta_v, rho_v;
    863:   Vec             f, X;
    864:   KSP             ksp;
    865:   PC              pc;
    866:   char            filename[PETSC_MAX_PATH_LEN];
    867:   DM              dms_quadrature, dms_mpoint;
    868:   PetscInt        nel, npe, npoints;
    869:   const PetscInt *element_list;
    870:   PetscInt        tk, nt, dump_freq;
    871:   PetscReal       dt, dt_max = 0.0;
    872:   PetscReal       vx[2], vy[2], max_v = 0.0, max_v_step, dh;
    873:   const char     *fieldnames[] = {"eta", "rho"};
    874:   Vec             pfields[2];
    875:   PetscInt        ppcell = 1;
    876:   PetscReal       time, delta_eta = 1.0;
    877:   PetscBool       randomize_coords = PETSC_FALSE;
    878:   PetscReal       randomize_fac    = 0.25;
    879:   PetscBool       no_view          = PETSC_FALSE;
    880:   PetscBool       isbddc;

    882:   PetscFunctionBeginUser;
    883:   /*
    884:     Generate the DMDA for the velocity and pressure spaces.
    885:     We use Q1 elements for both fields.
    886:     The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
    887:     The number of nodes in each direction is mx+1, my+1
    888:   */
    889:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
    890:   p_dof         = P_DOFS; /* p - pressure */
    891:   dof           = u_dof + p_dof;
    892:   stencil_width = 1;
    893:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &dm_stokes));
    894:   PetscCall(DMDASetElementType(dm_stokes, DMDA_ELEMENT_Q1));
    895:   PetscCall(DMSetMatType(dm_stokes, MATAIJ));
    896:   PetscCall(DMSetFromOptions(dm_stokes));
    897:   PetscCall(DMSetUp(dm_stokes));
    898:   PetscCall(DMDASetFieldName(dm_stokes, 0, "ux"));
    899:   PetscCall(DMDASetFieldName(dm_stokes, 1, "uy"));
    900:   PetscCall(DMDASetFieldName(dm_stokes, 2, "p"));

    902:   /* unit box [0,0.9142] x [0,1] */
    903:   PetscCall(DMDASetUniformCoordinates(dm_stokes, 0.0, 0.9142, 0.0, 1.0, 0., 0.));
    904:   dh = 1.0 / (PetscReal)mx;

    906:   /* Get local number of elements */
    907:   {
    908:     PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));

    910:     nel_local = nel;

    912:     PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
    913:   }

    915:   /* Create DMDA for representing scalar fields */
    916:   PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 1, &dm_coeff));

    918:   /* Create the swarm for storing quadrature point values */
    919:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_quadrature));
    920:   PetscCall(DMSetType(dms_quadrature, DMSWARM));
    921:   PetscCall(DMSetDimension(dms_quadrature, 2));
    922:   PetscCall(PetscObjectSetName((PetscObject)dms_quadrature, "Quadrature Swarm"));

    924:   /* Register fields for viscosity and density on the quadrature points */
    925:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "eta_q", 1, PETSC_REAL));
    926:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "rho_q", 1, PETSC_REAL));
    927:   PetscCall(DMSwarmFinalizeFieldRegister(dms_quadrature));
    928:   PetscCall(DMSwarmSetLocalSizes(dms_quadrature, nel_local * GAUSS_POINTS, 0));

    930:   /* Create the material point swarm */
    931:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_mpoint));
    932:   PetscCall(DMSetType(dms_mpoint, DMSWARM));
    933:   PetscCall(DMSetDimension(dms_mpoint, 2));
    934:   PetscCall(PetscObjectSetName((PetscObject)dms_mpoint, "Material Point Swarm"));

    936:   /* Configure the material point swarm to be of type Particle-In-Cell */
    937:   PetscCall(DMSwarmSetType(dms_mpoint, DMSWARM_PIC));

    939:   /*
    940:      Specify the DM to use for point location and projections
    941:      within the context of a PIC scheme
    942:   */
    943:   PetscCall(DMSwarmSetCellDM(dms_mpoint, dm_coeff));

    945:   /* Register fields for viscosity and density */
    946:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "eta", 1, PETSC_REAL));
    947:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "rho", 1, PETSC_REAL));
    948:   PetscCall(DMSwarmFinalizeFieldRegister(dms_mpoint));

    950:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL));
    951:   PetscCall(DMSwarmSetLocalSizes(dms_mpoint, nel_local * ppcell, 100));

    953:   /*
    954:     Layout the material points in space using the cell DM.
    955:     Particle coordinates are defined by cell wise using different methods.
    956:     - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
    957:                               corresponding to a Gauss quadrature rule with
    958:                               ppcell points in each direction.
    959:     - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
    960:                                 ppcell x ppcell quadralaterals defined within the
    961:                                 reference element.
    962:     - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
    963:                                     of each quadralateral obtained by sub-dividing
    964:                                     the reference element cell ppcell times.
    965:   */
    966:   PetscCall(DMSwarmInsertPointsUsingCellDM(dms_mpoint, DMSWARMPIC_LAYOUT_SUBDIVISION, ppcell));

    968:   /*
    969:     Defne a high resolution layer of material points across the material interface
    970:   */
    971:   {
    972:     PetscInt  npoints_dir_x[2];
    973:     PetscReal min[2], max[2];

    975:     npoints_dir_x[0] = (PetscInt)(0.9142 / (0.05 * dh));
    976:     npoints_dir_x[1] = (PetscInt)((0.25 - 0.15) / (0.05 * dh));
    977:     min[0]           = 0.0;
    978:     max[0]           = 0.9142;
    979:     min[1]           = 0.05;
    980:     max[1]           = 0.35;
    981:     PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
    982:   }

    984:   /*
    985:     Define a high resolution layer of material points near the surface of the domain
    986:     to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
    987:     when applied to buoyancy driven flow. The error in div(u) is O(h).
    988:   */
    989:   {
    990:     PetscInt  npoints_dir_x[2];
    991:     PetscReal min[2], max[2];

    993:     npoints_dir_x[0] = (PetscInt)(0.9142 / (0.25 * dh));
    994:     npoints_dir_x[1] = (PetscInt)(3.0 * dh / (0.25 * dh));
    995:     min[0]           = 0.0;
    996:     max[0]           = 0.9142;
    997:     min[1]           = 1.0 - 3.0 * dh;
    998:     max[1]           = 1.0 - 0.0001;
    999:     PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
    1000:   }

    1002:   PetscCall(DMView(dms_mpoint, PETSC_VIEWER_STDOUT_WORLD));

    1004:   /* Define initial material properties on each particle in the material point swarm */
    1005:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-delta_eta", &delta_eta, NULL));
    1006:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-randomize_coords", &randomize_coords, NULL));
    1007:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-randomize_fac", &randomize_fac, NULL));
    1008:   PetscCheck(randomize_fac <= 1.0, PETSC_COMM_WORLD, PETSC_ERR_USER, "The value of -randomize_fac should be <= 1.0");
    1009:   {
    1010:     PetscReal  *array_x, *array_e, *array_r;
    1011:     PetscInt    p;
    1012:     PetscRandom r;
    1013:     PetscMPIInt rank;

    1015:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

    1017:     PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
    1018:     PetscCall(PetscRandomSetInterval(r, -randomize_fac * dh, randomize_fac * dh));
    1019:     PetscCall(PetscRandomSetSeed(r, (unsigned long)rank));
    1020:     PetscCall(PetscRandomSeed(r));

    1022:     PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));

    1024:     /*
    1025:        Fetch the registered data from the material point DMSwarm.
    1026:        The fields "eta" and "rho" were registered by this example.
    1027:        The field identified by the variable DMSwarmPICField_coor
    1028:        was registered by the DMSwarm implementation when the function
    1029:          DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
    1030:        was called. The returned array defines the coordinates of each
    1031:        material point in the point swarm.
    1032:     */
    1033:     PetscCall(DMSwarmGetField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
    1034:     PetscCall(DMSwarmGetField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
    1035:     PetscCall(DMSwarmGetField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));

    1037:     PetscCall(DMSwarmGetLocalSize(dms_mpoint, &npoints));
    1038:     for (p = 0; p < npoints; p++) {
    1039:       PetscReal x_p[2], rr[2];

    1041:       if (randomize_coords) {
    1042:         PetscCall(PetscRandomGetValueReal(r, &rr[0]));
    1043:         PetscCall(PetscRandomGetValueReal(r, &rr[1]));
    1044:         array_x[2 * p + 0] += rr[0];
    1045:         array_x[2 * p + 1] += rr[1];
    1046:       }

    1048:       /* Get the coordinates of point, p */
    1049:       x_p[0] = array_x[2 * p + 0];
    1050:       x_p[1] = array_x[2 * p + 1];

    1052:       if (x_p[1] < (0.2 + 0.02 * PetscCosReal(PETSC_PI * x_p[0] / 0.9142))) {
    1053:         /* Material properties below the interface */
    1054:         array_e[p] = 1.0 * (1.0 / delta_eta);
    1055:         array_r[p] = 0.0;
    1056:       } else {
    1057:         /* Material properties above the interface */
    1058:         array_e[p] = 1.0;
    1059:         array_r[p] = 1.0;
    1060:       }
    1061:     }

    1063:     /*
    1064:        Restore the fetched data fields from the material point DMSwarm.
    1065:        Calling the Restore function invalidates the points array_r, array_e, array_x
    1066:        by setting them to NULL.
    1067:     */
    1068:     PetscCall(DMSwarmRestoreField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
    1069:     PetscCall(DMSwarmRestoreField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
    1070:     PetscCall(DMSwarmRestoreField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));

    1072:     PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
    1073:     PetscCall(PetscRandomDestroy(&r));
    1074:   }

    1076:   /*
    1077:      If the particle coordinates where randomly shifted, they may have crossed into another
    1078:      element, or into another sub-domain. To account for this we call the Migrate function.
    1079:   */
    1080:   if (randomize_coords) PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));

    1082:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
    1083:   if (!no_view) PetscCall(DMSwarmViewXDMF(dms_mpoint, "ic_coeff_dms.xmf"));

    1085:   /* project the swarm properties */
    1086:   PetscCall(DMCreateGlobalVector(dm_coeff, &pfields[0]));
    1087:   PetscCall(DMCreateGlobalVector(dm_coeff, &pfields[1]));
    1088:   PetscCall(DMSwarmProjectFields(dms_mpoint, NULL, 2, fieldnames, pfields, SCATTER_FORWARD));
    1089:   eta_v = pfields[0];
    1090:   rho_v = pfields[1];
    1091:   PetscCall(PetscObjectSetName((PetscObject)eta_v, "eta"));
    1092:   PetscCall(PetscObjectSetName((PetscObject)rho_v, "rho"));
    1093:   PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));

    1095:   /* view projected coefficients eta and rho */
    1096:   if (!no_view) {
    1097:     PetscViewer viewer;

    1099:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
    1100:     PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
    1101:     PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
    1102:     PetscCall(PetscViewerFileSetName(viewer, "ic_coeff_dmda.vts"));
    1103:     PetscCall(VecView(eta_v, viewer));
    1104:     PetscCall(VecView(rho_v, viewer));
    1105:     PetscCall(PetscViewerDestroy(&viewer));
    1106:   }

    1108:   PetscCall(DMCreateMatrix(dm_stokes, &A));
    1109:   PetscCall(DMCreateMatrix(dm_stokes, &B));
    1110:   PetscCall(DMCreateGlobalVector(dm_stokes, &f));
    1111:   PetscCall(DMCreateGlobalVector(dm_stokes, &X));

    1113:   PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
    1114:   PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
    1115:   PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));

    1117:   PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
    1118:   PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));

    1120:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
    1121:   PetscCall(KSPSetOptionsPrefix(ksp, "stokes_"));
    1122:   PetscCall(KSPSetDM(ksp, dm_stokes));
    1123:   PetscCall(KSPSetDMActive(ksp, KSP_DMACTIVE_ALL, PETSC_FALSE));
    1124:   PetscCall(KSPSetOperators(ksp, A, B));
    1125:   PetscCall(KSPSetFromOptions(ksp));
    1126:   PetscCall(KSPGetPC(ksp, &pc));
    1127:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &isbddc));
    1128:   if (isbddc) PetscCall(KSPSetOperators(ksp, A, A));

    1130:   /* Define u-v-p indices for fieldsplit */
    1131:   {
    1132:     PC             pc;
    1133:     const PetscInt ufields[] = {0, 1}, pfields[1] = {2};

    1135:     PetscCall(KSPGetPC(ksp, &pc));
    1136:     PetscCall(PCFieldSplitSetBlockSize(pc, 3));
    1137:     PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
    1138:     PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
    1139:   }

    1141:   /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
    1142:   {
    1143:     PC        pc, pc_u;
    1144:     KSP      *sub_ksp, ksp_u;
    1145:     PetscInt  nsplits;
    1146:     DM        dm_u;
    1147:     PetscBool is_pcfs;

    1149:     PetscCall(KSPGetPC(ksp, &pc));

    1151:     is_pcfs = PETSC_FALSE;
    1152:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs));

    1154:     if (is_pcfs) {
    1155:       PetscCall(KSPSetUp(ksp));
    1156:       PetscCall(KSPGetPC(ksp, &pc));
    1157:       PetscCall(PCFieldSplitGetSubKSP(pc, &nsplits, &sub_ksp));
    1158:       ksp_u = sub_ksp[0];
    1159:       PetscCall(PetscFree(sub_ksp));

    1161:       if (nsplits == 2) {
    1162:         PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 2, &dm_u));

    1164:         PetscCall(KSPSetDM(ksp_u, dm_u));
    1165:         PetscCall(KSPSetDMActive(ksp_u, KSP_DMACTIVE_ALL, PETSC_FALSE));
    1166:         PetscCall(DMDestroy(&dm_u));

    1168:         /* enforce galerkin coarse grids be used */
    1169:         PetscCall(KSPGetPC(ksp_u, &pc_u));
    1170:         PetscCall(PCMGSetGalerkin(pc_u, PC_MG_GALERKIN_PMAT));
    1171:       }
    1172:     }
    1173:   }

    1175:   dump_freq = 10;
    1176:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dump_freq", &dump_freq, NULL));
    1177:   nt = 10;
    1178:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nt", &nt, NULL));
    1179:   time = 0.0;
    1180:   for (tk = 1; tk <= nt; tk++) {
    1181:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... assemble\n"));
    1182:     PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
    1183:     PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
    1184:     PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));

    1186:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... bc imposition\n"));
    1187:     PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
    1188:     PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));

    1190:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... solve\n"));
    1191:     PetscCall(KSPSetOperators(ksp, A, isbddc ? A : B));
    1192:     PetscCall(KSPSolve(ksp, f, X));

    1194:     PetscCall(VecStrideMax(X, 0, NULL, &vx[1]));
    1195:     PetscCall(VecStrideMax(X, 1, NULL, &vy[1]));
    1196:     PetscCall(VecStrideMin(X, 0, NULL, &vx[0]));
    1197:     PetscCall(VecStrideMin(X, 1, NULL, &vy[0]));

    1199:     max_v_step = PetscMax(vx[0], vx[1]);
    1200:     max_v_step = PetscMax(max_v_step, vy[0]);
    1201:     max_v_step = PetscMax(max_v_step, vy[1]);
    1202:     max_v      = PetscMax(max_v, max_v_step);

    1204:     dt_max = 2.0;
    1205:     dt     = 0.5 * (dh / max_v_step);
    1206:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... max v %1.4e , dt %1.4e : [total] max v %1.4e , dt_max %1.4e\n", (double)max_v_step, (double)dt, (double)max_v, (double)dt_max));
    1207:     dt = PetscMin(dt_max, dt);

    1209:     /* advect */
    1210:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... advect\n"));
    1211:     PetscCall(MaterialPoint_AdvectRK1(dm_stokes, X, dt, dms_mpoint));

    1213:     /* migrate */
    1214:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... migrate\n"));
    1215:     PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));

    1217:     /* update cell population */
    1218:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... populate cells\n"));
    1219:     PetscCall(MaterialPoint_PopulateCell(dm_stokes, dms_mpoint));

    1221:     /* update coefficients on quadrature points */
    1222:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... project\n"));
    1223:     PetscCall(DMSwarmProjectFields(dms_mpoint, NULL, 2, fieldnames, pfields, SCATTER_FORWARD));
    1224:     eta_v = pfields[0];
    1225:     rho_v = pfields[1];
    1226:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... interp\n"));
    1227:     PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));

    1229:     if (tk % dump_freq == 0) {
    1230:       PetscViewer viewer;

    1232:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... write XDMF, VTS\n"));
    1233:       PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_coeff_dms.xmf", tk));
    1234:       PetscCall(DMSwarmViewXDMF(dms_mpoint, filename));

    1236:       PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_vp_dm.vts", tk));
    1237:       PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
    1238:       PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
    1239:       PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
    1240:       PetscCall(PetscViewerFileSetName(viewer, filename));
    1241:       PetscCall(VecView(X, viewer));
    1242:       PetscCall(PetscViewerDestroy(&viewer));
    1243:     }
    1244:     time += dt;
    1245:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step %" PetscInt_FMT " : time %1.2e \n", tk, (double)time));
    1246:   }

    1248:   PetscCall(KSPDestroy(&ksp));
    1249:   PetscCall(VecDestroy(&X));
    1250:   PetscCall(VecDestroy(&f));
    1251:   PetscCall(MatDestroy(&A));
    1252:   PetscCall(MatDestroy(&B));
    1253:   PetscCall(VecDestroy(&eta_v));
    1254:   PetscCall(VecDestroy(&rho_v));

    1256:   PetscCall(DMDestroy(&dms_mpoint));
    1257:   PetscCall(DMDestroy(&dms_quadrature));
    1258:   PetscCall(DMDestroy(&dm_coeff));
    1259:   PetscCall(DMDestroy(&dm_stokes));
    1260:   PetscFunctionReturn(PETSC_SUCCESS);
    1261: }

    1263: /*
    1264:  <sequential run>
    1265:  ./ex70 -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type lu -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type lu  -mx 80 -my 80  -stokes_ksp_converged_reason  -dump_freq 25  -stokes_ksp_rtol 1.0e-8 -build_twosided allreduce  -ppcell 2 -nt 4000 -delta_eta 1.0 -randomize_coords
    1266: */
    1267: int main(int argc, char **args)
    1268: {
    1269:   PetscInt  mx, my;
    1270:   PetscBool set = PETSC_FALSE;

    1272:   PetscFunctionBeginUser;
    1273:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
    1274:   mx = my = 10;
    1275:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
    1276:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
    1277:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mxy", &mx, &set));
    1278:   if (set) my = mx;
    1279:   PetscCall(SolveTimeDepStokes(mx, my));
    1280:   PetscCall(PetscFinalize());
    1281:   return 0;
    1282: }

    1284: /* -------------------------- helpers for boundary conditions -------------------------------- */
    1285: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
    1286: {
    1287:   DM                     cda;
    1288:   Vec                    coords;
    1289:   PetscInt               si, sj, nx, ny, i, j;
    1290:   PetscInt               M, N;
    1291:   DMDACoor2d           **_coords;
    1292:   const PetscInt        *g_idx;
    1293:   PetscInt              *bc_global_ids;
    1294:   PetscScalar           *bc_vals;
    1295:   PetscInt               nbcs;
    1296:   PetscInt               n_dofs;
    1297:   ISLocalToGlobalMapping ltogm;

    1299:   PetscFunctionBeginUser;
    1300:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1301:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1303:   PetscCall(DMGetCoordinateDM(da, &cda));
    1304:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1305:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1306:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1307:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1309:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
    1310:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

    1312:   /* init the entries to -1 so VecSetValues will ignore them */
    1313:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

    1315:   i = nx - 1;
    1316:   for (j = 0; j < ny; j++) {
    1317:     PetscInt local_id;

    1319:     local_id = i + j * nx;

    1321:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

    1323:     bc_vals[j] = 0.0;
    1324:   }
    1325:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1326:   nbcs = 0;
    1327:   if ((si + nx) == (M)) nbcs = ny;

    1329:   if (b) {
    1330:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1331:     PetscCall(VecAssemblyBegin(b));
    1332:     PetscCall(VecAssemblyEnd(b));
    1333:   }
    1334:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

    1336:   PetscCall(PetscFree(bc_vals));
    1337:   PetscCall(PetscFree(bc_global_ids));

    1339:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1340:   PetscFunctionReturn(PETSC_SUCCESS);
    1341: }

    1343: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
    1344: {
    1345:   DM                     cda;
    1346:   Vec                    coords;
    1347:   PetscInt               si, sj, nx, ny, i, j;
    1348:   PetscInt               M, N;
    1349:   DMDACoor2d           **_coords;
    1350:   const PetscInt        *g_idx;
    1351:   PetscInt              *bc_global_ids;
    1352:   PetscScalar           *bc_vals;
    1353:   PetscInt               nbcs;
    1354:   PetscInt               n_dofs;
    1355:   ISLocalToGlobalMapping ltogm;

    1357:   PetscFunctionBeginUser;
    1358:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1359:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1361:   PetscCall(DMGetCoordinateDM(da, &cda));
    1362:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1363:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1364:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1365:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1367:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
    1368:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

    1370:   /* init the entries to -1 so VecSetValues will ignore them */
    1371:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

    1373:   i = 0;
    1374:   for (j = 0; j < ny; j++) {
    1375:     PetscInt local_id;

    1377:     local_id = i + j * nx;

    1379:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

    1381:     bc_vals[j] = 0.0;
    1382:   }
    1383:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1384:   nbcs = 0;
    1385:   if (si == 0) nbcs = ny;

    1387:   if (b) {
    1388:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1389:     PetscCall(VecAssemblyBegin(b));
    1390:     PetscCall(VecAssemblyEnd(b));
    1391:   }

    1393:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

    1395:   PetscCall(PetscFree(bc_vals));
    1396:   PetscCall(PetscFree(bc_global_ids));

    1398:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1399:   PetscFunctionReturn(PETSC_SUCCESS);
    1400: }

    1402: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
    1403: {
    1404:   DM                     cda;
    1405:   Vec                    coords;
    1406:   PetscInt               si, sj, nx, ny, i, j;
    1407:   PetscInt               M, N;
    1408:   DMDACoor2d           **_coords;
    1409:   const PetscInt        *g_idx;
    1410:   PetscInt              *bc_global_ids;
    1411:   PetscScalar           *bc_vals;
    1412:   PetscInt               nbcs;
    1413:   PetscInt               n_dofs;
    1414:   ISLocalToGlobalMapping ltogm;

    1416:   PetscFunctionBeginUser;
    1417:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1418:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1420:   PetscCall(DMGetCoordinateDM(da, &cda));
    1421:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1422:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1423:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1424:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1426:   PetscCall(PetscMalloc1(nx, &bc_global_ids));
    1427:   PetscCall(PetscMalloc1(nx, &bc_vals));

    1429:   /* init the entries to -1 so VecSetValues will ignore them */
    1430:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

    1432:   j = ny - 1;
    1433:   for (i = 0; i < nx; i++) {
    1434:     PetscInt local_id;

    1436:     local_id = i + j * nx;

    1438:     bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];

    1440:     bc_vals[i] = 0.0;
    1441:   }
    1442:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1443:   nbcs = 0;
    1444:   if ((sj + ny) == (N)) nbcs = nx;

    1446:   if (b) {
    1447:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1448:     PetscCall(VecAssemblyBegin(b));
    1449:     PetscCall(VecAssemblyEnd(b));
    1450:   }
    1451:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, NULL, NULL));

    1453:   PetscCall(PetscFree(bc_vals));
    1454:   PetscCall(PetscFree(bc_global_ids));

    1456:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1457:   PetscFunctionReturn(PETSC_SUCCESS);
    1458: }

    1460: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
    1461: {
    1462:   DM                     cda;
    1463:   Vec                    coords;
    1464:   PetscInt               si, sj, nx, ny, i, j;
    1465:   PetscInt               M, N;
    1466:   DMDACoor2d           **_coords;
    1467:   const PetscInt        *g_idx;
    1468:   PetscInt              *bc_global_ids;
    1469:   PetscScalar           *bc_vals;
    1470:   PetscInt               nbcs;
    1471:   PetscInt               n_dofs;
    1472:   ISLocalToGlobalMapping ltogm;

    1474:   PetscFunctionBeginUser;
    1475:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1476:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1478:   PetscCall(DMGetCoordinateDM(da, &cda));
    1479:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1480:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1481:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1482:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1484:   PetscCall(PetscMalloc1(nx, &bc_global_ids));
    1485:   PetscCall(PetscMalloc1(nx, &bc_vals));

    1487:   /* init the entries to -1 so VecSetValues will ignore them */
    1488:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

    1490:   j = 0;
    1491:   for (i = 0; i < nx; i++) {
    1492:     PetscInt local_id;

    1494:     local_id = i + j * nx;

    1496:     bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];

    1498:     bc_vals[i] = 0.0;
    1499:   }
    1500:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1501:   nbcs = 0;
    1502:   if (sj == 0) nbcs = nx;

    1504:   if (b) {
    1505:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1506:     PetscCall(VecAssemblyBegin(b));
    1507:     PetscCall(VecAssemblyEnd(b));
    1508:   }
    1509:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

    1511:   PetscCall(PetscFree(bc_vals));
    1512:   PetscCall(PetscFree(bc_global_ids));

    1514:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1515:   PetscFunctionReturn(PETSC_SUCCESS);
    1516: }

    1518: /*
    1519:  Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
    1520:  Impose no slip boundray conditions on the top/bottom faces:   u_i n_i = 0, u_i t_i = 0
    1521: */
    1522: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes, Mat A, Vec f)
    1523: {
    1524:   PetscFunctionBeginUser;
    1525:   PetscCall(BCApplyZero_NORTH(dm_stokes, 0, A, f));
    1526:   PetscCall(BCApplyZero_NORTH(dm_stokes, 1, A, f));
    1527:   PetscCall(BCApplyZero_EAST(dm_stokes, 0, A, f));
    1528:   PetscCall(BCApplyZero_SOUTH(dm_stokes, 0, A, f));
    1529:   PetscCall(BCApplyZero_SOUTH(dm_stokes, 1, A, f));
    1530:   PetscCall(BCApplyZero_WEST(dm_stokes, 0, A, f));
    1531:   PetscFunctionReturn(PETSC_SUCCESS);
    1532: }

    1534: /*TEST

    1536:    test:
    1537:      suffix: 1
    1538:      args: -no_view
    1539:      requires: !complex double
    1540:      filter: grep -v atomic
    1541:      filter_output: grep -v atomic
    1542:    test:
    1543:      suffix: 1_matis
    1544:      requires: !complex double
    1545:      args: -no_view -dm_mat_type is
    1546:      filter: grep -v atomic
    1547:      filter_output: grep -v atomic
    1548:    testset:
    1549:      nsize: 4
    1550:      requires: !complex double
    1551:      args: -no_view -dm_mat_type is -stokes_ksp_type fetidp -mx 80 -my 80 -stokes_ksp_converged_reason -stokes_ksp_rtol 1.0e-8 -ppcell 2 -nt 4 -randomize_coords -stokes_ksp_error_if_not_converged
    1552:      filter: grep -v atomic
    1553:      filter_output: grep -v atomic
    1554:      test:
    1555:        suffix: fetidp
    1556:        args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
    1557:      test:
    1558:        suffix: fetidp_lumped
    1559:        args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd -stokes_fetidp_pc_lumped -stokes_fetidp_bddc_pc_bddc_dirichlet_pc_type none -stokes_fetidp_bddc_pc_bddc_switch_static
    1560:      test:
    1561:        suffix: fetidp_saddlepoint
    1562:        args: -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_norm_type natural -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1
    1563:      test:
    1564:        suffix: fetidp_saddlepoint_lumped
    1565:        args: -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_norm_type natural -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1 -stokes_fetidp_bddc_pc_bddc_dirichlet_pc_type none -stokes_fetidp_bddc_pc_bddc_switch_static -stokes_fetidp_pc_lumped
    1566: TEST*/