Actual source code: ex69.c

      1: #include <petscdt.h>
      2: #include <petscdraw.h>
      3: #include <petscviewer.h>
      4: #include <petscksp.h>
      5: #include <petscdmda.h>

      7: /*
      8:       Solves -Laplacian u = f,  u(-1) = u(1) = 0 with multiple spectral elements

     10:       Uses DMDA to manage the parallelization of the elements

     12:       This is not intended to be highly optimized in either memory usage or time, but strifes for simplicity.

     14: */

     16: typedef struct {
     17:   PetscInt   n;       /* number of nodes */
     18:   PetscReal *nodes;   /* GLL nodes */
     19:   PetscReal *weights; /* GLL weights */
     20: } PetscGLL;

     22: PetscErrorCode ComputeSolution(DM da, PetscGLL *gll, Vec u)
     23: {
     24:   PetscInt     j, xs, xn;
     25:   PetscScalar *uu, *xx;
     26:   PetscReal    xd;
     27:   Vec          x;

     29:   PetscFunctionBegin;
     30:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xn, NULL, NULL));
     31:   PetscCall(DMGetCoordinates(da, &x));
     32:   PetscCall(DMDAVecGetArray(da, x, &xx));
     33:   PetscCall(DMDAVecGetArray(da, u, &uu));
     34:   /* loop over local nodes */
     35:   for (j = xs; j < xs + xn; j++) {
     36:     xd    = xx[j];
     37:     uu[j] = (xd * xd - 1.0) * PetscCosReal(5. * PETSC_PI * xd);
     38:   }
     39:   PetscCall(DMDAVecRestoreArray(da, x, &xx));
     40:   PetscCall(DMDAVecRestoreArray(da, u, &uu));
     41:   PetscFunctionReturn(PETSC_SUCCESS);
     42: }

     44: /*
     45:       Evaluates \integral_{-1}^{1} f*v_i  where v_i is the ith basis polynomial via the GLL nodes and weights, since the v_i
     46:       basis function is zero at all nodes except the ith one the integral is simply the weight_i * f(node_i)
     47: */
     48: PetscErrorCode ComputeRhs(DM da, PetscGLL *gll, Vec b)
     49: {
     50:   PetscInt     i, j, xs, xn, n = gll->n;
     51:   PetscScalar *bb, *xx;
     52:   PetscReal    xd;
     53:   Vec          blocal, xlocal;

     55:   PetscFunctionBegin;
     56:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xn, NULL, NULL));
     57:   xs = xs / (n - 1);
     58:   xn = xn / (n - 1);
     59:   PetscCall(DMGetLocalVector(da, &blocal));
     60:   PetscCall(VecZeroEntries(blocal));
     61:   PetscCall(DMDAVecGetArray(da, blocal, &bb));
     62:   PetscCall(DMGetCoordinatesLocal(da, &xlocal));
     63:   PetscCall(DMDAVecGetArray(da, xlocal, &xx));
     64:   /* loop over local spectral elements */
     65:   for (j = xs; j < xs + xn; j++) {
     66:     /* loop over GLL points in each element */
     67:     for (i = 0; i < n; i++) {
     68:       xd = xx[j * (n - 1) + i];
     69:       bb[j * (n - 1) + i] += -gll->weights[i] * (-20. * PETSC_PI * xd * PetscSinReal(5. * PETSC_PI * xd) + (2. - (5. * PETSC_PI) * (5. * PETSC_PI) * (xd * xd - 1.)) * PetscCosReal(5. * PETSC_PI * xd));
     70:     }
     71:   }
     72:   PetscCall(DMDAVecRestoreArray(da, xlocal, &xx));
     73:   PetscCall(DMDAVecRestoreArray(da, blocal, &bb));
     74:   PetscCall(VecZeroEntries(b));
     75:   PetscCall(DMLocalToGlobalBegin(da, blocal, ADD_VALUES, b));
     76:   PetscCall(DMLocalToGlobalEnd(da, blocal, ADD_VALUES, b));
     77:   PetscCall(DMRestoreLocalVector(da, &blocal));
     78:   PetscFunctionReturn(PETSC_SUCCESS);
     79: }

     81: /*
     82:      Run with -build_twosided allreduce -pc_type bjacobi -sub_pc_type lu -q 16 -ksp_rtol 1.e-34 (or 1.e-14 for double precision)

     84:      -q <q> number of spectral elements to use
     85:      -N <N> maximum number of GLL points per element

     87: */
     88: int main(int argc, char **args)
     89: {
     90:   PetscGLL      gll;
     91:   PetscInt      N = 80, n, q = 8, xs, xn, j, l;
     92:   PetscReal   **A;
     93:   Mat           K;
     94:   KSP           ksp;
     95:   PC            pc;
     96:   Vec           x, b;
     97:   PetscInt     *rows;
     98:   PetscReal     norm, xc, yc, h;
     99:   PetscScalar  *f;
    100:   PetscDraw     draw;
    101:   PetscDrawLG   lg;
    102:   PetscDrawAxis axis;
    103:   DM            da;
    104:   PetscMPIInt   rank, size;

    106:   PetscFunctionBeginUser;
    107:   PetscCall(PetscInitialize(&argc, &args, NULL, NULL));
    108:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    109:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    110:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
    111:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-q", &q, NULL));

    113:   PetscCall(PetscDrawCreate(PETSC_COMM_WORLD, NULL, "Log(Error norm) vs Number of GLL points", 0, 0, 500, 500, &draw));
    114:   PetscCall(PetscDrawSetFromOptions(draw));
    115:   PetscCall(PetscDrawLGCreate(draw, 1, &lg));
    116:   PetscCall(PetscDrawLGSetUseMarkers(lg, PETSC_TRUE));
    117:   PetscCall(PetscDrawLGGetAxis(lg, &axis));
    118:   PetscCall(PetscDrawAxisSetLabels(axis, NULL, "Number of GLL points", "Log(Error Norm)"));

    120:   for (n = 4; n < N; n += 2) {
    121:     /*
    122:        da contains the information about the parallel layout of the elements
    123:     */
    124:     PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, q * (n - 1) + 1, 1, 1, NULL, &da));
    125:     PetscCall(DMSetFromOptions(da));
    126:     PetscCall(DMSetUp(da));
    127:     PetscCall(DMDAGetInfo(da, NULL, &q, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    128:     q = (q - 1) / (n - 1); /* number of spectral elements */

    130:     /*
    131:        gll simply contains the GLL node and weight values
    132:     */
    133:     PetscCall(PetscMalloc2(n, &gll.nodes, n, &gll.weights));
    134:     PetscCall(PetscDTGaussLobattoLegendreQuadrature(n, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, gll.nodes, gll.weights));
    135:     gll.n = n;
    136:     PetscCall(DMDASetGLLCoordinates(da, gll.n, gll.nodes));

    138:     /*
    139:        Creates the element stiffness matrix for the given gll
    140:     */
    141:     PetscCall(PetscGaussLobattoLegendreElementLaplacianCreate(gll.n, gll.nodes, gll.weights, &A));

    143:     /*
    144:       Scale the element stiffness and weights by the size of the element
    145:     */
    146:     h = 2.0 / q;
    147:     for (j = 0; j < n; j++) {
    148:       gll.weights[j] *= .5 * h;
    149:       for (l = 0; l < n; l++) A[j][l] = 2. * A[j][l] / h;
    150:     }

    152:     /*
    153:         Create the global stiffness matrix and add the element stiffness for each local element
    154:     */
    155:     PetscCall(DMCreateMatrix(da, &K));
    156:     PetscCall(MatSetOption(K, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
    157:     PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xn, NULL, NULL));
    158:     xs = xs / (n - 1);
    159:     xn = xn / (n - 1);
    160:     PetscCall(PetscMalloc1(n, &rows));
    161:     /*
    162:         loop over local elements
    163:     */
    164:     for (j = xs; j < xs + xn; j++) {
    165:       for (l = 0; l < n; l++) rows[l] = j * (n - 1) + l;
    166:       PetscCall(MatSetValues(K, n, rows, n, rows, &A[0][0], ADD_VALUES));
    167:     }
    168:     PetscCall(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY));
    169:     PetscCall(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY));

    171:     PetscCall(MatCreateVecs(K, &x, &b));
    172:     PetscCall(ComputeRhs(da, &gll, b));

    174:     /*
    175:         Replace the first and last rows/columns of the matrix with the identity to obtain the zero Dirichlet boundary conditions
    176:     */
    177:     rows[0] = 0;
    178:     rows[1] = q * (n - 1);
    179:     PetscCall(MatZeroRowsColumns(K, 2, rows, 1.0, x, b));
    180:     PetscCall(PetscFree(rows));

    182:     PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
    183:     PetscCall(KSPSetOperators(ksp, K, K));
    184:     PetscCall(KSPGetPC(ksp, &pc));
    185:     PetscCall(PCSetType(pc, PCLU));
    186:     PetscCall(KSPSetFromOptions(ksp));
    187:     PetscCall(KSPSolve(ksp, b, x));

    189:     /* compute the error to the continium problem */
    190:     PetscCall(ComputeSolution(da, &gll, b));
    191:     PetscCall(VecAXPY(x, -1.0, b));

    193:     /* compute the L^2 norm of the error */
    194:     PetscCall(VecGetArray(x, &f));
    195:     PetscCall(PetscGaussLobattoLegendreIntegrate(gll.n, gll.nodes, gll.weights, f, &norm));
    196:     PetscCall(VecRestoreArray(x, &f));
    197:     norm = PetscSqrtReal(norm);
    198:     PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "L^2 norm of the error %" PetscInt_FMT " %g\n", n, (double)norm));
    199:     PetscCheck(n <= 10 || norm <= 1.e-8, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Slower convergence than expected");
    200:     xc = (PetscReal)n;
    201:     yc = PetscLog10Real(norm);
    202:     PetscCall(PetscDrawLGAddPoint(lg, &xc, &yc));
    203:     PetscCall(PetscDrawLGDraw(lg));

    205:     PetscCall(VecDestroy(&b));
    206:     PetscCall(VecDestroy(&x));
    207:     PetscCall(KSPDestroy(&ksp));
    208:     PetscCall(MatDestroy(&K));
    209:     PetscCall(PetscGaussLobattoLegendreElementLaplacianDestroy(gll.n, gll.nodes, gll.weights, &A));
    210:     PetscCall(PetscFree2(gll.nodes, gll.weights));
    211:     PetscCall(DMDestroy(&da));
    212:   }
    213:   PetscCall(PetscDrawLGDestroy(&lg));
    214:   PetscCall(PetscDrawDestroy(&draw));
    215:   PetscCall(PetscFinalize());
    216:   return 0;
    217: }

    219: /*TEST

    221:   build:
    222:       requires: !complex

    224:    test:
    225:      requires: !single

    227:    test:
    228:      suffix: 2
    229:      nsize: 2
    230:      requires: superlu_dist

    232: TEST*/