Actual source code: ex68.c
petsc-3.13.6 2020-09-29
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*/