xref: /petsc/src/dm/tests/ex30.c (revision d71ae5a4db6382e7f06317b8d368875286fe9008)

static char help[] = "Tests DMSLICED operations\n\n";

#include <petscdmsliced.h>

int main(int argc, char *argv[])
{
  char         mat_type[256] = MATAIJ; /* default matrix type */
  MPI_Comm     comm;
  PetscMPIInt  rank, size;
  DM           slice;
  PetscInt     i, bs = 1, N = 5, n, m, rstart, ghosts[2], *d_nnz, *o_nnz, dfill[4] = {1, 0, 0, 1}, ofill[4] = {1, 1, 1, 1};
  PetscReal    alpha = 1, K = 1, rho0 = 1, u0 = 0, sigma = 0.2;
  PetscBool    useblock = PETSC_TRUE;
  PetscScalar *xx;
  Mat          A;
  Vec          x, b, lf;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, 0, help));
  comm = PETSC_COMM_WORLD;
  PetscCallMPI(MPI_Comm_size(comm, &size));
  PetscCallMPI(MPI_Comm_rank(comm, &rank));

  PetscOptionsBegin(comm, 0, "Options for DMSliced test", 0);
  {
    PetscCall(PetscOptionsInt("-n", "Global number of nodes", "", N, &N, NULL));
    PetscCall(PetscOptionsInt("-bs", "Block size (1 or 2)", "", bs, &bs, NULL));
    if (bs != 1) {
      PetscCheck(bs == 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Block size must be 1 or 2");
      PetscCall(PetscOptionsReal("-alpha", "Inverse time step for wave operator", "", alpha, &alpha, NULL));
      PetscCall(PetscOptionsReal("-K", "Bulk modulus of compressibility", "", K, &K, NULL));
      PetscCall(PetscOptionsReal("-rho0", "Reference density", "", rho0, &rho0, NULL));
      PetscCall(PetscOptionsReal("-u0", "Reference velocity", "", u0, &u0, NULL));
      PetscCall(PetscOptionsReal("-sigma", "Width of Gaussian density perturbation", "", sigma, &sigma, NULL));
      PetscCall(PetscOptionsBool("-block", "Use block matrix assembly", "", useblock, &useblock, NULL));
    }
    PetscCall(PetscOptionsString("-sliced_mat_type", "Matrix type to use (aij or baij)", "", mat_type, mat_type, sizeof(mat_type), NULL));
  }
  PetscOptionsEnd();

  /* Split ownership, set up periodic grid in 1D */
  n = PETSC_DECIDE;
  PetscCall(PetscSplitOwnership(comm, &n, &N));
  rstart = 0;
  PetscCallMPI(MPI_Scan(&n, &rstart, 1, MPIU_INT, MPI_SUM, comm));
  rstart -= n;
  ghosts[0] = (N + rstart - 1) % N;
  ghosts[1] = (rstart + n) % N;

  PetscCall(PetscMalloc2(n, &d_nnz, n, &o_nnz));
  for (i = 0; i < n; i++) {
    if (size > 1 && (i == 0 || i == n - 1)) {
      d_nnz[i] = 2;
      o_nnz[i] = 1;
    } else {
      d_nnz[i] = 3;
      o_nnz[i] = 0;
    }
  }
  PetscCall(DMSlicedCreate(comm, bs, n, 2, ghosts, d_nnz, o_nnz, &slice)); /* Currently does not copy X_nnz so we can't free them until after DMSlicedGetMatrix */

  if (!useblock) PetscCall(DMSlicedSetBlockFills(slice, dfill, ofill)); /* Irrelevant for baij formats */
  PetscCall(DMSetMatType(slice, mat_type));
  PetscCall(DMCreateMatrix(slice, &A));
  PetscCall(PetscFree2(d_nnz, o_nnz));
  PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));

  PetscCall(DMCreateGlobalVector(slice, &x));
  PetscCall(VecDuplicate(x, &b));

  PetscCall(VecGhostGetLocalForm(x, &lf));
  PetscCall(VecGetSize(lf, &m));
  PetscCheck(m == (n + 2) * bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size of local form %" PetscInt_FMT ", expected %" PetscInt_FMT, m, (n + 2) * bs);
  PetscCall(VecGetArray(lf, &xx));
  for (i = 0; i < n; i++) {
    PetscInt        row[2], col[9], im, ip;
    PetscScalar     v[12];
    const PetscReal xref = 2.0 * (rstart + i) / N - 1; /* [-1,1] */
    const PetscReal h    = 1.0 / N;                    /* grid spacing */
    im                   = (i == 0) ? n : i - 1;
    ip                   = (i == n - 1) ? n + 1 : i + 1;
    switch (bs) {
    case 1: /* Laplacian with periodic boundaries */
      col[0] = im;
      col[1] = i;
      col[2] = ip;
      v[0]   = -h;
      v[1]   = 2 * h;
      v[2]   = -h;
      PetscCall(MatSetValuesLocal(A, 1, &i, 3, col, v, INSERT_VALUES));
      xx[i] = PetscSinReal(xref * PETSC_PI);
      break;
    case 2: /* Linear acoustic wave operator in variables [rho, u], central differences, periodic, timestep 1/alpha */
      v[0]  = -0.5 * u0;
      v[1]  = -0.5 * K;
      v[2]  = alpha;
      v[3]  = 0;
      v[4]  = 0.5 * u0;
      v[5]  = 0.5 * K;
      v[6]  = -0.5 / rho0;
      v[7]  = -0.5 * u0;
      v[8]  = 0;
      v[9]  = alpha;
      v[10] = 0.5 / rho0;
      v[11] = 0.5 * u0;
      if (useblock) {
        row[0] = i;
        col[0] = im;
        col[1] = i;
        col[2] = ip;
        PetscCall(MatSetValuesBlockedLocal(A, 1, row, 3, col, v, INSERT_VALUES));
      } else {
        row[0] = 2 * i;
        row[1] = 2 * i + 1;
        col[0] = 2 * im;
        col[1] = 2 * im + 1;
        col[2] = 2 * i;
        col[3] = 2 * ip;
        col[4] = 2 * ip + 1;
        v[3]   = v[4];
        v[4]   = v[5]; /* pack values in first row */
        PetscCall(MatSetValuesLocal(A, 1, row, 5, col, v, INSERT_VALUES));
        col[2] = 2 * i + 1;
        v[8]   = v[9];
        v[9]   = v[10];
        v[10]  = v[11]; /* pack values in second row */
        PetscCall(MatSetValuesLocal(A, 1, row + 1, 5, col, v + 6, INSERT_VALUES));
      }
      /* Set current state (gaussian density perturbation) */
      xx[2 * i]     = 0.2 * PetscExpReal(-PetscSqr(xref) / (2 * PetscSqr(sigma)));
      xx[2 * i + 1] = 0;
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for block size %" PetscInt_FMT, bs);
    }
  }
  PetscCall(VecRestoreArray(lf, &xx));
  PetscCall(VecGhostRestoreLocalForm(x, &lf));
  PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

  PetscCall(MatMult(A, x, b));
  PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(VecView(b, PETSC_VIEWER_STDOUT_WORLD));

  /* Update the ghosted values, view the result on rank 0. */
  PetscCall(VecGhostUpdateBegin(b, INSERT_VALUES, SCATTER_FORWARD));
  PetscCall(VecGhostUpdateEnd(b, INSERT_VALUES, SCATTER_FORWARD));
  if (rank == 0) {
    PetscCall(VecGhostGetLocalForm(b, &lf));
    PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Local form of b on rank 0, last two nodes are ghost nodes\n"));
    PetscCall(VecView(lf, PETSC_VIEWER_STDOUT_SELF));
    PetscCall(VecGhostRestoreLocalForm(b, &lf));
  }

  PetscCall(DMDestroy(&slice));
  PetscCall(VecDestroy(&x));
  PetscCall(VecDestroy(&b));
  PetscCall(MatDestroy(&A));
  PetscCall(PetscFinalize());
  return 0;
}

/*TEST

   test:
      nsize: 2
      args: -bs 2 -block 0 -sliced_mat_type baij -alpha 10 -u0 0.1

   test:
      suffix: 2
      nsize: 2
      args: -bs 2 -block 1 -sliced_mat_type aij -alpha 10 -u0 0.1

   test:
      suffix: 3
      nsize: 2
      args: -bs 2 -block 0 -sliced_mat_type aij -alpha 10 -u0 0.1

TEST*/