[petsc-dev] Bug in PCSetUp_FieldSplit/PCReset_FieldSplit
Sander Arens
Sander.Arens at ugent.be
Tue Nov 22 03:10:38 CST 2016
It works now for the particular case when there's a SNESReset before a
SNESSolve, but when you put a SNESReset and SNESSolve again after the first
SNESSolve it breaks again (I noticed this first playing around with
TSArkimex when it broke at this line
<https://bitbucket.org/petsc/petsc/src/7538c81a3e5819d80e753a9b491da681486bdef7/src/ts/impls/arkimex/arkimex.c?at=master&fileviewer=file-view-default#arkimex.c-761>
).
Destroying the fieldsplit linked list in PCReset_FieldSplit seemed to work
(see diff file)...
On 21 November 2016 at 21:02, Matthew Knepley <knepley at gmail.com> wrote:
> Okay, try out the fix in next.
>
> Thanks,
>
> Matt
>
> On Thu, Nov 17, 2016 at 3:11 PM, Matthew Knepley <knepley at gmail.com>
> wrote:
>
>> On Thu, Nov 17, 2016 at 1:22 PM, Sander Arens <Sander.Arens at ugent.be>
>> wrote:
>>
>>> I think there's a bug in PCSetUp_FieldSplit or PCReset_FieldSplit. I can
>>> get this by putting a SNESReset before the SNESSolve in snes ex77 and
>>> running ./config/builder2.py check src/snes/examples/tutorials/ex77.c.
>>> Can somebody reproduce this?
>>>
>>
>> You are right. Fixing...
>>
>> Matt
>>
>>
>>> Thanks,
>>> Sander
>>>
>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20161122/d91fe061/attachment.html>
-------------- next part --------------
diff --git a/master:src/ksp/pc/impls/fieldsplit/fieldsplit.c b/sanderarens/fix-fieldsplit-reset:src/ksp/pc/impls/fieldsplit/fieldsplit.c
index 9ae7608..97f263d 100644
--- a/master:src/ksp/pc/impls/fieldsplit/fieldsplit.c
+++ b/sanderarens/fix-fieldsplit-reset:src/ksp/pc/impls/fieldsplit/fieldsplit.c
@@ -16,7 +16,7 @@ struct _PC_FieldSplitLink {
PetscInt nfields;
PetscInt *fields,*fields_col;
VecScatter sctx;
- IS is,is_col,is_orig;
+ IS is,is_col;
PC_FieldSplitLink next,previous;
PetscLogEvent event;
};
@@ -44,7 +44,6 @@ typedef struct {
KSP kspschur; /* The solver for S */
KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
PC_FieldSplitLink head;
- PetscBool reset; /* indicates PCReset() has been last called on this object, hack */
PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
@@ -324,10 +323,10 @@ static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
PetscFunctionBegin;
/*
- Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset.
+ Kinda messy, but at least this now uses DMCreateFieldDecomposition().
Should probably be rewritten.
*/
- if (!ilink || jac->reset) {
+ if (!ilink) {
ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);CHKERRQ(ierr);
if (pc->dm && jac->dm_splits && !stokes && !coupling) {
@@ -407,14 +406,8 @@ static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr);
ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr);
ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr);
- if (jac->reset) {
- if (!jac->head) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"At reset jac must have head from previous setup");
- jac->head->is = rest;
- jac->head->next->is = zerodiags;
- } else {
- ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr);
- ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr);
- }
+ ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr);
+ ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr);
ierr = ISDestroy(&zerodiags);CHKERRQ(ierr);
ierr = ISDestroy(&rest);CHKERRQ(ierr);
} else if (coupling) {
@@ -425,18 +418,11 @@ static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
ierr = MatFindOffBlockDiagonalEntries(pc->mat,&coupling);CHKERRQ(ierr);
ierr = ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);CHKERRQ(ierr);
ierr = ISSetIdentity(rest);CHKERRQ(ierr);
- if (jac->reset) {
- if (!jac->head) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"At reset jac must have head from previous setup");
- jac->head->is = rest;
- jac->head->next->is = coupling;
- } else {
- ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr);
- ierr = PCFieldSplitSetIS(pc,"1",coupling);CHKERRQ(ierr);
- }
+ ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr);
+ ierr = PCFieldSplitSetIS(pc,"1",coupling);CHKERRQ(ierr);
ierr = ISDestroy(&coupling);CHKERRQ(ierr);
ierr = ISDestroy(&rest);CHKERRQ(ierr);
} else {
- if (jac->reset && !jac->isrestrict) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used");
ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr);
if (!fieldsplit_default) {
/* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
@@ -1189,20 +1175,23 @@ static PetscErrorCode PCReset_FieldSplit(PC pc)
PetscFunctionBegin;
while (ilink) {
+ ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr);
ierr = KSPReset(ilink->ksp);CHKERRQ(ierr);
ierr = VecDestroy(&ilink->x);CHKERRQ(ierr);
ierr = VecDestroy(&ilink->y);CHKERRQ(ierr);
ierr = VecDestroy(&ilink->z);CHKERRQ(ierr);
ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr);
- if (!ilink->is_orig) { /* save the original IS */
- ierr = PetscObjectReference((PetscObject)ilink->is);CHKERRQ(ierr);
- ilink->is_orig = ilink->is;
- }
ierr = ISDestroy(&ilink->is);CHKERRQ(ierr);
ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr);
+ ierr = PetscFree(ilink->splitname);CHKERRQ(ierr);
+ ierr = PetscFree(ilink->fields);CHKERRQ(ierr);
+ ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr);
next = ilink->next;
ilink = next;
+ ierr = PetscFree(ilink);CHKERRQ(ierr);
}
+ jac->head = NULL;
+ jac->nsplits = 0;
ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr);
if (jac->mat && jac->mat != jac->pmat) {
ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr);
@@ -1220,7 +1209,6 @@ static PetscErrorCode PCReset_FieldSplit(PC pc)
ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr);
ierr = MatDestroy(&jac->B);CHKERRQ(ierr);
ierr = MatDestroy(&jac->C);CHKERRQ(ierr);
- jac->reset = PETSC_TRUE;
jac->isrestrict = PETSC_FALSE;
PetscFunctionReturn(0);
}
@@ -1235,17 +1223,6 @@ static PetscErrorCode PCDestroy_FieldSplit(PC pc)
PetscFunctionBegin;
ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr);
- while (ilink) {
- ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr);
- ierr = ISDestroy(&ilink->is_orig);CHKERRQ(ierr);
- next = ilink->next;
- ierr = PetscFree(ilink->splitname);CHKERRQ(ierr);
- ierr = PetscFree(ilink->fields);CHKERRQ(ierr);
- ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr);
- ierr = PetscFree(ilink);CHKERRQ(ierr);
- ilink = next;
- }
- ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr);
ierr = PetscFree(pc->data);CHKERRQ(ierr);
ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr);
ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr);
@@ -1447,11 +1424,7 @@ static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
while(ilink) {
IS isrl,isr;
PC subpc;
- if (jac->reset) {
- ierr = ISEmbed(ilink->is_orig, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr);
- } else {
- ierr = ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr);
- }
+ ierr = ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr);
ierr = ISGetLocalSize(isrl,&localsize);CHKERRQ(ierr);
ierr = PetscMalloc1(localsize,&indc);CHKERRQ(ierr);
ierr = ISGetIndices(isrl,&ind);CHKERRQ(ierr);
@@ -1472,15 +1445,9 @@ static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
if(flg) {
IS iszl,isz;
MPI_Comm comm;
- if (jac->reset) {
- ierr = ISGetLocalSize(ilink->is_orig,&localsize);CHKERRQ(ierr);
- comm = PetscObjectComm((PetscObject)ilink->is_orig);
- ierr = ISEmbed(isy, ilink->is_orig, PETSC_TRUE, &iszl);CHKERRQ(ierr);
- } else {
- ierr = ISGetLocalSize(ilink->is,&localsize);CHKERRQ(ierr);
- comm = PetscObjectComm((PetscObject)ilink->is);
- ierr = ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);CHKERRQ(ierr);
- }
+ ierr = ISGetLocalSize(ilink->is,&localsize);CHKERRQ(ierr);
+ comm = PetscObjectComm((PetscObject)ilink->is);
+ ierr = ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);CHKERRQ(ierr);
ierr = MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr);
sizez -= localsize;
ierr = ISGetLocalSize(iszl,&localsize);CHKERRQ(ierr);
More information about the petsc-dev
mailing list