Actual source code: ex68.c

  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;

 23:   PetscFunctionBeginUser;
 24:   PetscCall(MatMult(A, x, f));
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

 28: PetscErrorCode ComputeJacobianLinear(SNES snes, Vec x, Mat A, Mat J, void *ctx)
 29: {
 30:   PetscFunctionBeginUser;
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

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

 38:   PetscFunctionBeginUser;
 39:   PetscCall(VecSet(b, -3.0));
 40:   PetscCall(MatGetOwnershipRange(A, &rStart, &rEnd));
 41:   for (row = rStart; row < rEnd; ++row) {
 42:     PetscScalar val = 1.0;

 44:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES));
 45:   }
 46:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 47:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 48:   PetscFunctionReturn(PETSC_SUCCESS);
 49: }

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

 56:   PetscFunctionBeginUser;
 57:   PetscCall(VecDuplicate(b, &errorVec));
 58:   PetscCall(VecWAXPY(errorVec, -1.0, b, u));
 59:   PetscCall(VecNorm(errorVec, NORM_2, &error));
 60:   PetscCall(VecNorm(b, NORM_2, &norm));
 61:   PetscCheck(error / norm <= 1000. * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm));
 62:   PetscCall(VecDestroy(&errorVec));
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: PetscErrorCode ConstructProblem2(Mat A, Vec b)
 67: {
 68:   PetscInt N = 10, constraintSize = 4;
 69:   PetscInt row;

 71:   PetscFunctionBeginUser;
 72:   PetscCall(VecSet(b, -3.0));
 73:   for (row = 0; row < constraintSize; ++row) {
 74:     PetscScalar vals[2] = {1.0, 1.0};
 75:     PetscInt    cols[2];

 77:     cols[0] = row;
 78:     cols[1] = row + N - constraintSize;
 79:     PetscCall(MatSetValues(A, 1, &row, 2, cols, vals, INSERT_VALUES));
 80:   }
 81:   for (row = constraintSize; row < N - constraintSize; ++row) {
 82:     PetscScalar val = 1.0;

 84:     PetscCall(MatSetValues(A, 1, &row, 1, &row, &val, INSERT_VALUES));
 85:   }
 86:   for (row = N - constraintSize; row < N; ++row) {
 87:     PetscInt    col = row - (N - constraintSize);
 88:     PetscScalar val = 1.0;

 90:     PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES));
 91:   }
 92:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 93:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: PetscErrorCode CheckProblem2(Mat A, Vec b, Vec u)
 98: {
 99:   PetscInt           N = 10, constraintSize = 4, r;
100:   PetscReal          norm, error;
101:   const PetscScalar *uArray, *bArray;

103:   PetscFunctionBeginUser;
104:   PetscCall(VecNorm(b, NORM_2, &norm));
105:   PetscCall(VecGetArrayRead(u, &uArray));
106:   PetscCall(VecGetArrayRead(b, &bArray));
107:   error = 0.0;
108:   for (r = 0; r < constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r + N - constraintSize]));

110:   PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm));
111:   error = 0.0;
112:   for (r = constraintSize; r < N - constraintSize; ++r) error += PetscRealPart(PetscSqr(uArray[r] - bArray[r]));

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

118:   PetscCheck(error / norm <= 10000 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Relative error %g is too large", (double)(error / norm));
119:   PetscCall(VecRestoreArrayRead(u, &uArray));
120:   PetscCall(VecRestoreArrayRead(b, &bArray));
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }

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

132:   PetscFunctionBeginUser;
133:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
134:   comm = PETSC_COMM_WORLD;
135:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-problem", &problem, NULL));
136:   PetscCall(VecCreate(comm, &u));
137:   PetscCall(VecSetSizes(u, PETSC_DETERMINE, N));
138:   PetscCall(VecSetFromOptions(u));
139:   PetscCall(VecDuplicate(u, &r));
140:   PetscCall(VecDuplicate(u, &b));

142:   PetscCall(MatCreate(comm, &A));
143:   PetscCall(MatSetSizes(A, PETSC_DETERMINE, PETSC_DETERMINE, N, N));
144:   PetscCall(MatSetFromOptions(A));
145:   PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
146:   J = A;

148:   switch (problem) {
149:   case 1:
150:     PetscCall(ConstructProblem1(A, b));
151:     break;
152:   case 2:
153:     PetscCall(ConstructProblem2(A, b));
154:     break;
155:   default:
156:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem);
157:   }

159:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
160:   PetscCall(SNESSetJacobian(snes, A, J, ComputeJacobianLinear, NULL));
161:   PetscCall(SNESSetFunction(snes, r, ComputeFunctionLinear, A));
162:   PetscCall(SNESSetFromOptions(snes));

164:   PetscCall(SNESSolve(snes, b, u));
165:   PetscCall(VecView(u, NULL));

167:   switch (problem) {
168:   case 1:
169:     PetscCall(CheckProblem1(A, b, u));
170:     break;
171:   case 2:
172:     PetscCall(CheckProblem2(A, b, u));
173:     break;
174:   default:
175:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid problem number %" PetscInt_FMT, problem);
176:   }

178:   if (A != J) PetscCall(MatDestroy(&A));
179:   PetscCall(MatDestroy(&J));
180:   PetscCall(VecDestroy(&u));
181:   PetscCall(VecDestroy(&r));
182:   PetscCall(VecDestroy(&b));
183:   PetscCall(SNESDestroy(&snes));
184:   PetscCall(PetscFinalize());
185:   return 0;
186: }

188: /*TEST

190:    test:
191:      args: -snes_monitor

193:    test:
194:      suffix: 2
195:      args: -problem 2 -pc_type jacobi -snes_monitor

197: TEST*/