Actual source code: ex68.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      1: static char help[] = "Test problems for Schur complement solvers.\n\n\n";

      3:  #include <petscsnes.h>

      5: /*
      6: Test 1:
      7:   I u = b

      9:   solution: u = b

     11: Test 2:
     12:   / I 0 I \  / u_1 \   / b_1 \
     13:   | 0 I 0 | |  u_2 | = | b_2 |
     14:   \ I 0 0 /  \ u_3 /   \ b_3 /

     16:   solution: u_1 = b_3, u_2 = b_2, u_3 = b_1 - b_3
     17: */

     19: PetscErrorCode ComputeFunctionLinear(SNES snes, Vec x, Vec f, void *ctx)
     20: {
     21:   Mat            A = (Mat) ctx;

     25:   MatMult(A, x, f);
     26:   return(0);
     27: }

     29: PetscErrorCode ComputeJacobianLinear(SNES snes, Vec x, Mat A, Mat J, void *ctx)
     30: {
     32:   return(0);
     33: }

     35: PetscErrorCode ConstructProblem1(Mat A, Vec b)
     36: {
     37:   PetscInt       rStart, rEnd, row;

     41:   VecSet(b, -3.0);
     42:   MatGetOwnershipRange(A, &rStart, &rEnd);
     43:   for (row = rStart; row < rEnd; ++row) {
     44:     PetscScalar val = 1.0;

     46:     MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES);
     47:   }
     48:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
     49:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
     50:   return(0);
     51: }

     53: PetscErrorCode CheckProblem1(Mat A, Vec b, Vec u)
     54: {
     55:   Vec            errorVec;
     56:   PetscReal      norm, error;

     60:   VecDuplicate(b, &errorVec);
     61:   VecWAXPY(errorVec, -1.0, b, u);
     62:   VecNorm(errorVec, NORM_2, &error);
     63:   VecNorm(b, NORM_2, &norm);
     64:   if (error/norm > 1000.*PETSC_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
     65:   VecDestroy(&errorVec);
     66:   return(0);
     67: }

     69: PetscErrorCode ConstructProblem2(Mat A, Vec b)
     70: {
     71:   PetscInt       N = 10, constraintSize = 4;
     72:   PetscInt       row;

     76:   VecSet(b, -3.0);
     77:   for (row = 0; row < constraintSize; ++row) {
     78:     PetscScalar vals[2] = {1.0, 1.0};
     79:     PetscInt    cols[2];

     81:     cols[0] = row; cols[1] = row + N - constraintSize;
     82:     MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES);
     83:   }
     84:   for (row = constraintSize; row < N - constraintSize; ++row) {
     85:     PetscScalar val = 1.0;

     87:     MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES);
     88:   }
     89:   for (row = N - constraintSize; row < N; ++row) {
     90:     PetscInt    col = row - (N - constraintSize);
     91:     PetscScalar val = 1.0;

     93:     MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES);
     94:   }
     95:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
     96:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
     97:   return(0);
     98: }

    100: PetscErrorCode CheckProblem2(Mat A, Vec b, Vec u)
    101: {
    102:   PetscInt          N = 10, constraintSize = 4, r;
    103:   PetscReal         norm, error;
    104:   const PetscScalar *uArray, *bArray;
    105:   PetscErrorCode    ierr;

    108:   VecNorm(b, NORM_2, &norm);
    109:   VecGetArrayRead(u, &uArray);
    110:   VecGetArrayRead(b, &bArray);
    111:   error = 0.0;
    112:   for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N-constraintSize]));

    114:   if (error/norm > 10000*PETSC_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
    115:   error = 0.0;
    116:   for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r]));

    118:   if (error/norm > 10000*PETSC_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
    119:   error = 0.0;
    120:   for (r = N - constraintSize; r < N; ++r) error += PetscRealPart(PetscSqr(uArray[r] - (bArray[r - (N-constraintSize)] - bArray[r])));

    122:   if (error/norm > 10000*PETSC_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", error/norm);
    123:   VecRestoreArrayRead(u, &uArray);
    124:   VecRestoreArrayRead(b, &bArray);
    125:   return(0);
    126: }

    128: int main(int argc, char **argv)
    129: {
    130:   MPI_Comm       comm;
    131:   SNES           snes;                 /* nonlinear solver */
    132:   Vec            u,r,b;                /* solution, residual, and rhs vectors */
    133:   Mat            A,J;                  /* Jacobian matrix */
    134:   PetscInt       problem = 1, N = 10;

    137:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
    138:   comm = PETSC_COMM_WORLD;
    139:   PetscOptionsGetInt(NULL,NULL, "-problem", &problem, NULL);
    140:   VecCreate(comm, &u);
    141:   VecSetSizes(u, PETSC_DETERMINE, N);
    142:   VecSetFromOptions(u);
    143:   VecDuplicate(u, &r);
    144:   VecDuplicate(u, &b);

    146:   MatCreate(comm, &A);
    147:   MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N);
    148:   MatSetFromOptions(A);
    149:   MatSeqAIJSetPreallocation(A, 5, NULL);
    150:   J    = A;

    152:   switch (problem) {
    153:   case 1:
    154:     ConstructProblem1(A, b);
    155:     break;
    156:   case 2:
    157:     ConstructProblem2(A, b);
    158:     break;
    159:   default:
    160:     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem);
    161:   }

    163:   SNESCreate(PETSC_COMM_WORLD, &snes);
    164:   SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL);
    165:   SNESSetFunction(snes, r, ComputeFunctionLinear, A);
    166:   SNESSetFromOptions(snes);

    168:   SNESSolve(snes, b, u);
    169:   VecView(u, NULL);

    171:   switch (problem) {
    172:   case 1:
    173:     CheckProblem1(A, b, u);
    174:     break;
    175:   case 2:
    176:     CheckProblem2(A, b, u);
    177:     break;
    178:   default:
    179:     SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %d", problem);
    180:   }

    182:   if (A != J) {
    183:     MatDestroy(&A);
    184:   }
    185:   MatDestroy(&J);
    186:   VecDestroy(&u);
    187:   VecDestroy(&r);
    188:   VecDestroy(&b);
    189:   SNESDestroy(&snes);
    190:   PetscFinalize();
    191:   return ierr;
    192: }

    195: /*TEST

    197:    test:
    198:      args: -snes_monitor

    200:    test:
    201:      suffix: 2
    202:      args: -problem 2 -pc_type jacobi -snes_monitor

    204: TEST*/