Lines Matching refs:s
97 PetscErrorCode StokesSetupPC(Stokes *s, KSP ksp) in StokesSetupPC() argument
105 PetscCall(PCFieldSplitSetIS(pc, "0", s->isg[0])); in StokesSetupPC()
106 PetscCall(PCFieldSplitSetIS(pc, "1", s->isg[1])); in StokesSetupPC()
107 if (s->userPC) PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS)); in StokesSetupPC()
108 if (s->userKSP) { in StokesSetupPC()
111 PetscCall(KSPSetOperators(subksp[1], s->myS, s->myS)); in StokesSetupPC()
117 PetscErrorCode StokesWriteSolution(Stokes *s) in StokesWriteSolution() argument
128 PetscCall(VecGetArrayRead(s->x, &array)); in StokesWriteSolution()
131 for (j = 0; j < s->ny; j++) { in StokesWriteSolution()
132 for (i = 0; i < s->nx; i++) { in StokesWriteSolution()
133 n = j * s->nx + i; in StokesWriteSolution()
134 …s->hx + s->hx / 2), (double)(j * s->hy + s->hy / 2), (double)PetscRealPart(array[n]), (double)Pets… in StokesWriteSolution()
137 PetscCall(VecRestoreArrayRead(s->x, &array)); in StokesWriteSolution()
143 PetscErrorCode StokesSetupIndexSets(Stokes *s) in StokesSetupIndexSets() argument
147 PetscCall(MatNestGetISs(s->A, s->isg, NULL)); in StokesSetupIndexSets()
151 PetscErrorCode StokesSetupVectors(Stokes *s) in StokesSetupVectors() argument
155 PetscCall(VecCreate(PETSC_COMM_WORLD, &s->x)); in StokesSetupVectors()
156 PetscCall(VecSetSizes(s->x, PETSC_DECIDE, 3 * s->nx * s->ny)); in StokesSetupVectors()
157 PetscCall(VecSetType(s->x, VECMPI)); in StokesSetupVectors()
160 PetscCall(VecDuplicate(s->x, &s->y)); in StokesSetupVectors()
161 PetscCall(StokesExactSolution(s)); in StokesSetupVectors()
164 PetscCall(VecDuplicate(s->x, &s->b)); in StokesSetupVectors()
165 PetscCall(StokesRhs(s)); in StokesSetupVectors()
169 PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j) in StokesGetPosition() argument
175 n = row % (s->nx * s->ny); in StokesGetPosition()
176 *i = n % s->nx; in StokesGetPosition()
177 *j = (n - (*i)) / s->nx; in StokesGetPosition()
181 PetscErrorCode StokesExactSolution(Stokes *s) in StokesExactSolution() argument
189 PetscCall(VecGetSubVector(s->y, s->isg[0], &y0)); in StokesExactSolution()
192 PetscCall(StokesGetPosition(s, row, &i, &j)); in StokesExactSolution()
193 if (row < s->nx * s->ny) { in StokesExactSolution()
194 val = StokesExactVelocityX(j * s->hy + s->hy / 2); in StokesExactSolution()
202 PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0)); in StokesExactSolution()
205 PetscCall(VecGetSubVector(s->y, s->isg[1], &y1)); in StokesExactSolution()
208 PetscCall(StokesGetPosition(s, row, &i, &j)); in StokesExactSolution()
209 val = StokesExactPressure(i * s->hx + s->hx / 2); in StokesExactSolution()
214 PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1)); in StokesExactSolution()
218 PetscErrorCode StokesRhs(Stokes *s) in StokesRhs() argument
226 PetscCall(VecGetSubVector(s->b, s->isg[0], &b0)); in StokesRhs()
229 PetscCall(StokesGetPosition(s, row, &i, &j)); in StokesRhs()
230 if (row < s->nx * s->ny) { in StokesRhs()
231 PetscCall(StokesRhsMomX(s, i, j, &val)); in StokesRhs()
233 PetscCall(StokesRhsMomY(s, i, j, &val)); in StokesRhs()
237 PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0)); in StokesRhs()
240 PetscCall(VecGetSubVector(s->b, s->isg[1], &b1)); in StokesRhs()
243 PetscCall(StokesGetPosition(s, row, &i, &j)); in StokesRhs()
244 PetscCall(StokesRhsMass(s, i, j, &val)); in StokesRhs()
245 if (s->matsymmetric) val = -1.0 * val; in StokesRhs()
248 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1)); in StokesRhs()
252 PetscErrorCode StokesSetupMatBlock00(Stokes *s) in StokesSetupMatBlock00() argument
260 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[0])); in StokesSetupMatBlock00()
261 PetscCall(MatSetOptionsPrefix(s->subA[0], "a00_")); in StokesSetupMatBlock00()
262 …PetscCall(MatSetSizes(s->subA[0], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, 2 * s->nx * s->ny… in StokesSetupMatBlock00()
263 PetscCall(MatSetType(s->subA[0], MATMPIAIJ)); in StokesSetupMatBlock00()
264 PetscCall(MatMPIAIJSetPreallocation(s->subA[0], 5, NULL, 5, NULL)); in StokesSetupMatBlock00()
265 PetscCall(MatGetOwnershipRange(s->subA[0], &start, &end)); in StokesSetupMatBlock00()
268 PetscCall(StokesGetPosition(s, row, &i, &j)); in StokesSetupMatBlock00()
270 PetscCall(StokesStencilLaplacian(s, i, j, &sz, cols, vals)); in StokesSetupMatBlock00()
272 if (row >= s->nx * s->ny) { in StokesSetupMatBlock00()
273 for (i = 0; i < sz; i++) cols[i] += s->nx * s->ny; in StokesSetupMatBlock00()
276 PetscCall(MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES)); in StokesSetupMatBlock00()
278 PetscCall(MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY)); in StokesSetupMatBlock00()
279 PetscCall(MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY)); in StokesSetupMatBlock00()
283 PetscErrorCode StokesSetupMatBlock01(Stokes *s) in StokesSetupMatBlock01() argument
291 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[1])); in StokesSetupMatBlock01()
292 PetscCall(MatSetOptionsPrefix(s->subA[1], "a01_")); in StokesSetupMatBlock01()
293 PetscCall(MatSetSizes(s->subA[1], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, s->nx * s->ny)); in StokesSetupMatBlock01()
294 PetscCall(MatSetType(s->subA[1], MATMPIAIJ)); in StokesSetupMatBlock01()
295 PetscCall(MatMPIAIJSetPreallocation(s->subA[1], 5, NULL, 5, NULL)); in StokesSetupMatBlock01()
296 PetscCall(MatGetOwnershipRange(s->subA[1], &start, &end)); in StokesSetupMatBlock01()
298 PetscCall(MatSetOption(s->subA[1], MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); in StokesSetupMatBlock01()
301 PetscCall(StokesGetPosition(s, row, &i, &j)); in StokesSetupMatBlock01()
303 if (row < s->nx * s->ny) { in StokesSetupMatBlock01()
304 PetscCall(StokesStencilGradientX(s, i, j, &sz, cols, vals)); in StokesSetupMatBlock01()
306 PetscCall(StokesStencilGradientY(s, i, j, &sz, cols, vals)); in StokesSetupMatBlock01()
308 PetscCall(MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES)); in StokesSetupMatBlock01()
310 PetscCall(MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY)); in StokesSetupMatBlock01()
311 PetscCall(MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY)); in StokesSetupMatBlock01()
315 PetscErrorCode StokesSetupMatBlock10(Stokes *s) in StokesSetupMatBlock10() argument
319 PetscCall(MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2])); in StokesSetupMatBlock10()
320 if (!s->matsymmetric) PetscCall(MatScale(s->subA[2], -1.0)); in StokesSetupMatBlock10()
321 PetscCall(MatSetOptionsPrefix(s->subA[2], "a10_")); in StokesSetupMatBlock10()
325 PetscErrorCode StokesSetupMatBlock11(Stokes *s) in StokesSetupMatBlock11() argument
329 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[3])); in StokesSetupMatBlock11()
330 PetscCall(MatSetOptionsPrefix(s->subA[3], "a11_")); in StokesSetupMatBlock11()
331 PetscCall(MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx * s->ny, s->nx * s->ny)); in StokesSetupMatBlock11()
332 PetscCall(MatSetType(s->subA[3], MATMPIAIJ)); in StokesSetupMatBlock11()
333 PetscCall(MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL)); in StokesSetupMatBlock11()
334 PetscCall(MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY)); in StokesSetupMatBlock11()
335 PetscCall(MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY)); in StokesSetupMatBlock11()
339 PetscErrorCode StokesSetupApproxSchur(Stokes *s) in StokesSetupApproxSchur() argument
351 PetscCall(VecSetSizes(diag, PETSC_DECIDE, 2 * s->nx * s->ny)); in StokesSetupApproxSchur()
353 PetscCall(MatGetDiagonal(s->subA[0], diag)); in StokesSetupApproxSchur()
357 PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); /* (*warning* overwrites subA[1]) */ in StokesSetupApproxSchur()
358 PetscCall(MatMatMult(s->subA[2], s->subA[1], MAT_INITIAL_MATRIX, PETSC_DETERMINE, &s->myS)); in StokesSetupApproxSchur()
359 PetscCall(MatScale(s->myS, -1.0)); in StokesSetupApproxSchur()
362 PetscCall(MatGetDiagonal(s->subA[0], diag)); in StokesSetupApproxSchur()
363 PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); in StokesSetupApproxSchur()
368 PetscErrorCode StokesSetupMatrix(Stokes *s) in StokesSetupMatrix() argument
371 PetscCall(StokesSetupMatBlock00(s)); in StokesSetupMatrix()
372 PetscCall(StokesSetupMatBlock01(s)); in StokesSetupMatrix()
373 PetscCall(StokesSetupMatBlock10(s)); in StokesSetupMatrix()
374 PetscCall(StokesSetupMatBlock11(s)); in StokesSetupMatrix()
375 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A)); in StokesSetupMatrix()
376 PetscCall(StokesSetupApproxSchur(s)); in StokesSetupMatrix()
380 PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *co… in StokesStencilLaplacian() argument
382 PetscInt p = j * s->nx + i, w = p - 1, e = p + 1, s2 = p - s->nx, n = p + s->nx; in StokesStencilLaplacian()
383 PetscScalar ae = s->hy / s->hx, aeb = 0; in StokesStencilLaplacian()
384 PetscScalar aw = s->hy / s->hx, awb = s->hy / (s->hx / 2); in StokesStencilLaplacian()
385 PetscScalar as = s->hx / s->hy, asb = s->hx / (s->hy / 2); in StokesStencilLaplacian()
386 PetscScalar an = s->hx / s->hy, anb = s->hx / (s->hy / 2); in StokesStencilLaplacian()
397 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */ in StokesStencilLaplacian()
405 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */ in StokesStencilLaplacian()
413 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */ in StokesStencilLaplacian()
431 } else if (i == s->nx - 1) { /* east boundary */ in StokesStencilLaplacian()
451 } else if (j == s->ny - 1) { /* north boundary */ in StokesStencilLaplacian()
477 PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *co… in StokesStencilGradientX() argument
479 PetscInt p = j * s->nx + i, w = p - 1, e = p + 1; in StokesStencilGradientX()
480 PetscScalar ae = s->hy / 2, aeb = s->hy; in StokesStencilGradientX()
481 PetscScalar aw = -s->hy / 2, awb = 0; in StokesStencilGradientX()
490 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */ in StokesStencilGradientX()
496 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */ in StokesStencilGradientX()
502 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */ in StokesStencilGradientX()
514 } else if (i == s->nx - 1) { /* east boundary */ in StokesStencilGradientX()
528 } else if (j == s->ny - 1) { /* north boundary */ in StokesStencilGradientX()
548 PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *co… in StokesStencilGradientY() argument
550 PetscInt p = j * s->nx + i, s2 = p - s->nx, n = p + s->nx; in StokesStencilGradientY()
551 PetscScalar as = -s->hx / 2, asb = 0; in StokesStencilGradientY()
552 PetscScalar an = s->hx / 2, anb = 0; in StokesStencilGradientY()
561 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */ in StokesStencilGradientY()
567 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */ in StokesStencilGradientY()
573 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */ in StokesStencilGradientY()
587 } else if (i == s->nx - 1) { /* east boundary */ in StokesStencilGradientY()
601 } else if (j == s->ny - 1) { /* north boundary */ in StokesStencilGradientY()
619 PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) in StokesRhsMomX() argument
621 PetscScalar y = j * s->hy + s->hy / 2; in StokesRhsMomX()
622 PetscScalar awb = s->hy / (s->hx / 2); in StokesRhsMomX()
633 PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) in StokesRhsMomY() argument
640 PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) in StokesRhsMass() argument
642 PetscScalar y = j * s->hy + s->hy / 2; in StokesRhsMass()
643 PetscScalar aeb = s->hy; in StokesRhsMass()
654 PetscErrorCode StokesCalcResidual(Stokes *s) in StokesCalcResidual() argument
661 PetscCall(VecScale(s->b, -1.0)); in StokesCalcResidual()
662 PetscCall(MatMultAdd(s->A, s->x, s->b, s->b)); in StokesCalcResidual()
665 PetscCall(VecGetSubVector(s->b, s->isg[0], &b0)); in StokesCalcResidual()
668 PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0)); in StokesCalcResidual()
671 PetscCall(VecGetSubVector(s->b, s->isg[1], &b1)); in StokesCalcResidual()
674 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1)); in StokesCalcResidual()
677 PetscCall(VecNorm(s->b, NORM_2, &val)); in StokesCalcResidual()
682 PetscErrorCode StokesCalcError(Stokes *s) in StokesCalcError() argument
684 PetscScalar scale = PetscSqrtReal((double)s->nx * s->ny); in StokesCalcError()
690 PetscCall(VecAXPY(s->y, -1.0, s->x)); in StokesCalcError()
693 PetscCall(VecGetSubVector(s->y, s->isg[0], &y0)); in StokesCalcError()
696 PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0)); in StokesCalcError()
699 PetscCall(VecGetSubVector(s->y, s->isg[1], &y1)); in StokesCalcError()
702 PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1)); in StokesCalcError()
705 PetscCall(VecNorm(s->y, NORM_2, &val)); in StokesCalcError()
712 Stokes s; in main() local
717 s.nx = 4; in main()
718 s.ny = 6; in main()
719 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &s.nx, NULL)); in main()
720 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ny", &s.ny, NULL)); in main()
721 s.hx = 2.0 / s.nx; in main()
722 s.hy = 1.0 / s.ny; in main()
723 s.matsymmetric = PETSC_FALSE; in main()
724 PetscCall(PetscOptionsGetBool(NULL, NULL, "-mat_set_symmetric", &s.matsymmetric, NULL)); in main()
725 s.userPC = s.userKSP = PETSC_FALSE; in main()
726 PetscCall(PetscOptionsHasName(NULL, NULL, "-user_pc", &s.userPC)); in main()
727 PetscCall(PetscOptionsHasName(NULL, NULL, "-user_ksp", &s.userKSP)); in main()
729 PetscCall(StokesSetupMatrix(&s)); in main()
730 PetscCall(StokesSetupIndexSets(&s)); in main()
731 PetscCall(StokesSetupVectors(&s)); in main()
734 PetscCall(KSPSetOperators(ksp, s.A, s.A)); in main()
736 PetscCall(StokesSetupPC(&s, ksp)); in main()
737 PetscCall(KSPSolve(ksp, s.b, s.x)); in main()
740 PetscCall(StokesCalcResidual(&s)); in main()
741 PetscCall(StokesCalcError(&s)); in main()
742 PetscCall(StokesWriteSolution(&s)); in main()
745 PetscCall(MatDestroy(&s.subA[0])); in main()
746 PetscCall(MatDestroy(&s.subA[1])); in main()
747 PetscCall(MatDestroy(&s.subA[2])); in main()
748 PetscCall(MatDestroy(&s.subA[3])); in main()
749 PetscCall(MatDestroy(&s.A)); in main()
750 PetscCall(VecDestroy(&s.x)); in main()
751 PetscCall(VecDestroy(&s.b)); in main()
752 PetscCall(VecDestroy(&s.y)); in main()
753 PetscCall(MatDestroy(&s.myS)); in main()