[petsc-users] MatSetValuesStencil for cols not in the dmda stencil width

[Matthew Knepley]

On Thu, Jan 26, 2017 at 10:20 AM, Xiangdong <epscodes at gmail.com> wrote:

> Thanks, Matt.  For a cell in DMDA 3d with global id (ix,iy,iz), what is
> the global row id of that cell corresponding to the matrix generated by
> DMCreateMatrix? It is not always ix + iy*Nx + iz*Ny*Nx for
> different da_processors_xyz, is it? How can I obtain that global row id?
>

It is not simple:

 1) Figure out the process division using


http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMDA/DMDAGetInfo.html

 2) Figure out which process p your i,j,k in on

  You have to march through the processes, knowing how big the chunk of
indices is on each,
  which you get from 1) and


http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMDA/DMDAGetOwnershipRanges.html

 3) Add up the sizes from all process less than p, and then the offset on
process p. This the global offset

This is not a scalable thing to do, which is why we do not include it in
the API.

   Matt


> Thanks.
>
> Xiangdong
>
>
>
> On Wed, Jan 25, 2017 at 10:53 AM, Matthew Knepley <knepley at gmail.com>
> wrote:
>
>> On Wed, Jan 25, 2017 at 9:17 AM, Xiangdong <epscodes at gmail.com> wrote:
>>
>>> Hello everyone,
>>>
>>> I have a question on setting matrix entries which are not in the stencil
>>> width. Take ksp ex45.c as an example,
>>> http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examp
>>> les/tutorials/ex45.c.html
>>>
>>> Instead of setting up the standard 7-point stencil, now for each cell,
>>> the matrix also has a additional dependency on the cell (Mx, My, Mz).
>>> Namely, for each row, the col corresponding to (Mx, My, Mz) is always
>>> nonzero. I modify the example code to add this entries:
>>>
>>> +  MatSetOption(B,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
>>> +  MatSetOption(jac,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
>>>
>>> +  v[7] = 100; col[7].i = mx-1; col[7].j = my-1; col[7].k = mz-1;
>>> +  ierr = MatSetValuesStencil(B,1,&row,8,col,v,INSERT_VALUES);CHKERRQ(
>>> ierr);
>>>
>>> It is okay to for np=1, but crash for np>=2 with the error message:
>>>
>>
>> You can do this, but
>>
>> 1) You cannot use MatSetStencil, since your entry is not actually in your
>> stencil. You will have to make a separate call to MatSetValues() using the
>> global index.
>>
>> 2) The nonzero pattern we have allocated will be wrong. You will have to
>> set the MatOption which gives an error on new nonzeros to PETSC_FALSE.
>>
>> 3) You will have a dense row in your Jacobian, which is hard to make
>> perform well, and also stymies most preconditioners.
>>
>>   Thanks,
>>
>>     Matt
>>
>>
>>> [0]PETSC ERROR: --------------------- Error Message
>>> --------------------------------------------------------------
>>> [0]PETSC ERROR: Argument out of range
>>> [0]PETSC ERROR: Local index 342 too large 244 (max) at 7
>>>
>>> [0]PETSC ERROR: #1 ISLocalToGlobalMappingApply() line 423 in
>>> petsc-3.7.4/src/vec/is/utils/isltog.c
>>> [0]PETSC ERROR: #2 MatSetValuesLocal() line 2052 in
>>> petsc-3.7.4/src/mat/interface/matrix.c
>>> [0]PETSC ERROR: #3 MatSetValuesStencil() line 1447 in
>>> petsc-3.7.4/src/mat/interface/matrix.c
>>> [0]PETSC ERROR: #4 ComputeMatrix() line 151 in extest.c
>>>
>>> Can I add new entries to the cols not in the stencil width into the dmda
>>> matrix or Jacobian?
>>>
>>> Attached please find the modifed ex45 example, the diff file as well as
>>> the run log.
>>>
>>> Thanks for your help.
>>>
>>> Best,
>>> Xiangdong
>>>
>>>
>>
>>
>> --
>> 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-users/attachments/20170126/4a59b9e5/attachment.html>