Actual source code: init.cxx

      1: #include <petscdmda.h>
      2: #include <adolc/adalloc.h>

      4: /*
      5:    REQUIRES configuration of PETSc with option --download-adolc.

      7:    For documentation on ADOL-C, see
      8:      $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
      9: */

     11: /*
     12:   Wrapper function for allocating contiguous memory in a 2d array

     14:   Input parameters:
     15:   m,n - number of rows and columns of array, respectively

     17:   Output parameter:
     18:   A   - pointer to array for which memory is allocated

     20:   Note: Only arrays of doubles are currently accounted for in ADOL-C's myalloc2 function.
     21: */
     22: template <class T>
     23: PetscErrorCode AdolcMalloc2(PetscInt m, PetscInt n, T **A[])
     24: {
     25:   PetscFunctionBegin;
     26:   *A = myalloc2(m, n);
     27:   PetscFunctionReturn(PETSC_SUCCESS);
     28: }

     30: /*
     31:   Wrapper function for freeing memory allocated with AdolcMalloc2

     33:   Input parameter:
     34:   A - array to free memory of

     36:   Note: Only arrays of doubles are currently accounted for in ADOL-C's myfree2 function.
     37: */
     38: template <class T>
     39: PetscErrorCode AdolcFree2(T **A)
     40: {
     41:   PetscFunctionBegin;
     42:   myfree2(A);
     43:   PetscFunctionReturn(PETSC_SUCCESS);
     44: }

     46: /*
     47:   Shift indices in an array of type T to endow it with ghost points.
     48:   (e.g. This works for arrays of adoubles or arrays (of structs) thereof.)

     50:   Input parameters:
     51:   da   - distributed array upon which variables are defined
     52:   cgs  - contiguously allocated 1-array with as many entries as there are
     53:          interior and ghost points, in total

     55:   Output parameter:
     56:   array - contiguously allocated array of the appropriate dimension with
     57:           ghost points, pointing to the 1-array
     58: */
     59: template <class T>
     60: PetscErrorCode GiveGhostPoints(DM da, T *cgs, void *array)
     61: {
     62:   PetscInt dim;

     64:   PetscFunctionBegin;
     65:   PetscCall(DMDAGetInfo(da, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
     66:   if (dim == 1) {
     67:     PetscCall(GiveGhostPoints1d(da, (T **)array));
     68:   } else if (dim == 2) {
     69:     PetscCall(GiveGhostPoints2d(da, cgs, (T ***)array));
     70:   } else PetscCheck(dim != 3, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "GiveGhostPoints3d not yet implemented"); // TODO
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: /*
     75:   Shift indices in a 1-array of type T to endow it with ghost points.
     76:   (e.g. This works for arrays of adoubles or arrays (of structs) thereof.)

     78:   Input parameters:
     79:   da  - distributed array upon which variables are defined

     81:   Output parameter:
     82:   a1d - contiguously allocated 1-array
     83: */
     84: template <class T>
     85: PetscErrorCode GiveGhostPoints1d(DM da, T *a1d[])
     86: {
     87:   PetscInt gxs;

     89:   PetscFunctionBegin;
     90:   PetscCall(DMDAGetGhostCorners(da, &gxs, NULL, NULL, NULL, NULL, NULL));
     91:   *a1d -= gxs;
     92:   PetscFunctionReturn(PETSC_SUCCESS);
     93: }

     95: /*
     96:   Shift indices in a 2-array of type T to endow it with ghost points.
     97:   (e.g. This works for arrays of adoubles or arrays (of structs) thereof.)

     99:   Input parameters:
    100:   da  - distributed array upon which variables are defined
    101:   cgs - contiguously allocated 1-array with as many entries as there are
    102:         interior and ghost points, in total

    104:   Output parameter:
    105:   a2d - contiguously allocated 2-array with ghost points, pointing to the
    106:         1-array
    107: */
    108: template <class T>
    109: PetscErrorCode GiveGhostPoints2d(DM da, T *cgs, T **a2d[])
    110: {
    111:   PetscInt gxs, gys, gxm, gym;

    113:   PetscFunctionBegin;
    114:   PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gxm, &gym, NULL));
    115:   for (PetscInt j = 0; j < gym; j++) (*a2d)[j] = cgs + j * gxm - gxs;
    116:   *a2d -= gys;
    117:   PetscFunctionReturn(PETSC_SUCCESS);
    118: }

    120: /*
    121:   Create a rectangular sub-identity of the m x m identity matrix, as an array.

    123:   Input parameters:
    124:   n - number of (adjacent) rows to take in slice
    125:   s - starting row index

    127:   Output parameter:
    128:   S - resulting n x m submatrix
    129: */
    130: template <class T>
    131: PetscErrorCode Subidentity(PetscInt n, PetscInt s, T **S)
    132: {
    133:   PetscInt i;

    135:   PetscFunctionBegin;
    136:   for (i = 0; i < n; i++) S[i][i + s] = 1.;
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: /*
    141:   Create an identity matrix, as an array.

    143:   Input parameter:
    144:   n - number of rows/columns
    145:   I - n x n array with memory pre-allocated
    146: */
    147: template <class T>
    148: PetscErrorCode Identity(PetscInt n, T **I)
    149: {
    150:   PetscFunctionBegin;
    151:   PetscCall(Subidentity(n, 0, I));
    152:   PetscFunctionReturn(PETSC_SUCCESS);
    153: }