[petsc-users] fieldsplit: Unhandled case, must have at least two fields, not 1!

[Nico Schlömer]

Hi Jed,

thanks for the quick answer!
I can of course generate the strided index sets (01x34x6...,
xx2xx5x...) myself and insert them via PCFieldSplitSetIS(), but I
thought that ultimately does the same thing that I would have liked
PCFieldSplitSetFields() to do for me. Not true?

--Nico

On Thu, Nov 28, 2013 at 4:25 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
> Nico Schlömer <nico.schloemer at gmail.com> writes:
>
>> Hi all,
>>
>> I'm toying around with PETSc's field split capabilities for the
>> 2D-Stokes problem, and similar to ex42 I call
>>
>>   // set of the operators
>>   // set of a bunch of options
>>   PetscOptionsSetValue("-pc_type", "fieldsplit");
>>   PC pc;
>>   ierr = KSPGetPC(*_ksp, &pc);
>>   const PetscInt ufields[] = {0,1}, pfields[] = {2};
>>   ierr = PCFieldSplitSetBlockSize(pc,3);
>>   PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
>>   PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
>>
>> When finally solving, however, I'm getting the error message
>>
>> [0]PETSC ERROR: Petsc has generated inconsistent data!
>> [0]PETSC ERROR: Unhandled case, must have at least two fields, not 1!
>
> Short answer: use PCFieldSplitSetIS().  From the man page:
>
>     Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field.
>
>      The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
>      size is three then one can define a field as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean
>      0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
>      where the numbered entries indicate what is in the field.
>
>      This function is called once per split (it creates a new split each time).  Solve options
>      for this split will be available under the prefix -fieldsplit_SPLITNAME_.
>
>      Developer Note: This routine does not actually create the IS representing the split, that is delayed
>      until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be
>      available when this routine is called.