fieldsplit.c (ab57588789c9da5248d756019688f1b680cba327) fieldsplit.c (21635b76ed28834bc274bfdfd6051f93dfbfcb4f)
1
2#include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
3#include <petscdm.h>
4
5/*
6 There is a nice discussion of block preconditioners in
7
8[El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier–Stokes equations

--- 606 unchanged lines hidden (view full) ---

615 if (kspA != kspInner) {
616 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr);
617 }
618 if (kspUpper != kspA) {
619 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr);
620 }
621 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr);
622 } else {
1
2#include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
3#include <petscdm.h>
4
5/*
6 There is a nice discussion of block preconditioners in
7
8[El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier–Stokes equations

--- 606 unchanged lines hidden (view full) ---

615 if (kspA != kspInner) {
616 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr);
617 }
618 if (kspUpper != kspA) {
619 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr);
620 }
621 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr);
622 } else {
623 KSP ksp;
624 const char *Dprefix;
625 char schurprefix[256];
626 char schurtestoption[256];
627 MatNullSpace sp;
628 PetscBool flg;
629
630 /* extract the A01 and A10 matrices */
631 ilink = jac->head;
632 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
633 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr);
634 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
635 ilink = ilink->next;
636 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
637 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr);
638 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
639
640 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
641 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr);
642 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr);
623 const char *Dprefix;
624 char schurprefix[256];
625 char schurtestoption[256];
626 MatNullSpace sp;
627 PetscBool flg;
628
629 /* extract the A01 and A10 matrices */
630 ilink = jac->head;
631 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
632 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr);
633 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
634 ilink = ilink->next;
635 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr);
636 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr);
637 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
638
639 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
640 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr);
641 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr);
643 ierr = MatSchurComplementGetKSP(jac->schur, &ksp);CHKERRQ(ierr);
644 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr);
645 /* Indent this deeper to emphasize the "inner" nature of this solver. */
646 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject) pc, 2);CHKERRQ(ierr);
647 ierr = KSPSetOptionsPrefix(ksp, schurprefix);CHKERRQ(ierr);
648 ierr = MatSchurComplementSet(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr);
649
650 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr);
651 if (sp) {
652 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);
653 }
654
655 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr);
656 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr);
657 if (flg) {
642 ierr = MatSchurComplementSet(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr);
643
644 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr);
645 if (sp) {
646 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);
647 }
648
649 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr);
650 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr);
651 if (flg) {
658 DM dmInner;
652 DM dmInner;
653 KSP kspInner;
659
654
655 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr);
656 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr);
657 /* Indent this deeper to emphasize the "inner" nature of this solver. */
658 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr);
659 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr);
660
660 /* Set DM for new solver */
661 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr);
661 /* Set DM for new solver */
662 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr);
662 ierr = KSPSetDM(ksp, dmInner);CHKERRQ(ierr);
663 ierr = KSPSetDMActive(ksp, PETSC_FALSE);CHKERRQ(ierr);
664 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr);
665 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr);
663 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr);
664 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr);
666 } else {
665 } else {
667 ierr = MatSchurComplementSetKSP(jac->schur, jac->head->ksp);CHKERRQ(ierr);
666 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
667 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
668 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
669 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
670 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
671 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
672 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr);
673 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr);
668 }
674 }
675 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr);
676 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr);
669 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr);
670
671 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr);
672 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr);
673 if (flg) {
674 DM dmInner;
675
676 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr);

--- 1186 unchanged lines hidden --- 		677 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr);
678
679 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr);
680 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr);
681 if (flg) {
682 DM dmInner;
683
684 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr);

--- 1186 unchanged lines hidden ---