[petsc-users] Field wise residual norm

[Dharmendar Reddy]

Hello,
        Is the following code on the lines of the approach you suggested:

How do i handle the fortran call to DMCreateFieldIS where i am not
interested  in fieldNames and numFields as output.
Do i pass PETSC_NULL_OBJECT ?

subroutine computeResidualNormV2(meshDM,F,fieldId,fNorm)
  implicit none
#include "finclude/petsc.h90"
  ! IO Variables
  DM, intent(inout) :: meshDM
  Vec, intent(inout) :: F
  integer, intent(in ) :: fieldId
  real(WP), intent(inout) :: fNorm

  ! Local Variables
  PetscErrorCode :: ierr
  integer :: numField
  IS :: subfieldIS
  DM :: subDM
  character(len=50), pointer :: fieldName(:)
  Is, pointer :: fieldIS(:)
  VecScatter :: vscat

  call DMCreateSubDM(meshDM,1, fieldId, subfieldIS,subdm, ierr)

  call DMCreateFieldIS(meshDM, numField, fieldName, fieldIS, ierr)

  call DMGetGlobalVector(subDM, gVecsubField, ierr)

  call VecScatterCreate(F,fieldIS(fieldId),gVecsubField,subfieldIS,vscat,ierr)

  call VecScatterBegin(vscat, F, gVecsubField, INSET_VALUES,
SCATTER_FORWARD, ierr)
  call VecScatterEnd(vscat, F, gVecsubField, INSET_VALUES,
SCATTER_FORWARD, ierr)

  call VecNorm( gVecsubField, NORM_2, fNorm, ierr)

  call DMRestoreGlobalVector(subDM, gVecsubField, ierr)

end subroutine computeResidualNormV2

On Mon, Jan 13, 2014 at 4:53 PM, Dharmendar Reddy
<dharmareddy84 at gmail.com> wrote:
> I may need to read more about those functions to see how i can get the
> norms. But a quick look seem to suggest me that i will need lesser
> function calls than what i am thinking of doing right now.
>
> This is how my approach...Am i doing it right ?
>
> subroutine computeResidualNorm(meshDM,F,fieldId,fNorm)
>   implicit none
> #include "finclude/petsc.h90"
>   ! IO Variables
>   DM, intent(inout) :: meshDM
>   Vec, intent(inout) :: F
>   integer, intent(in ) :: fieldId
>   real(WP), intent(inout) :: fNorm
>
>   ! Local Variables
>   PetscErrorCode :: ierr
>   PetscSection :: Section
>   PetscScalar, pointer :: FArray(:)
>   Vec :: gVec, lVec
>   PetscScalar, pointer :: lArray(:)
>   integer :: numFieldDof,offset,dofId
>   integer :: CellId,CellIdStart,CellIdEnd
>   PetscScalar :: zero
>
>   zero = 0.0_WP
>   call DMGetGlobalVector(meshDM, gVec, ierr)
>   call DMGetLocalVector(meshDM, lVec, ierr)
>   call Vecset(lVec, zero, ierr)
>
>   ! Get acess to default section
>   call DMGetDefaultSection(meshDM,section,ierr)
>   ! Get chart
>   call  PetscSectionGetChart(section,cellIdstart,CellIdEnd,ierr)
>   call vecGetArrayF90(F, FArray, ierr)
>   call VecGetArrayF90(lVec, lArray, ierr)
>   do cellId=cellIdstart,cellIdEnd-1
>     call PetscSectionGetFieldDof(section,CellId,fieldId,numFieldDof,ierr)
>     if(numFieldDof > 0) then
>       call PetscSectionGetFieldOffset(section,cellId,fieldId,offset,ierr)
>       do dofId=0,numFieldDof-1
>         ic = offset+dofId
>         lArray(ic) = FArray(ic)
>       end do
>     end if
>   end do
>   call VecRestoreArrayF90(lVec, lArray, ierr)
>   call vecRestoreArrayF90(F, Farray,ierr)
>
>   call DMLocalToGlobalBegin(meshDM, lVec, INSERT_VALUES, gVec, ierr)
>   call DMLocalToGlobalEnd(meshDM, lVec, INSERT_VALUES, gVec, ierr)
>
>   call DMRestoreLocalVector(meshDM, lVec, ierr)
>
>   call VecNorm( gVec, NORM_2, fNorm, ierr)
>
>   call DMRestoreGlobalVector(meshDM, gVec, ierr)
> end subroutine computeResidualNorm
>
> On Mon, Jan 13, 2014 at 4:15 PM, Matthew Knepley <knepley at gmail.com> wrote:
>> On Mon, Jan 13, 2014 at 4:11 PM, Dharmendar Reddy <dharmareddy84 at gmail.com>
>> wrote:
>>>
>>> Hello,
>>>         Is there a petsc function which can help accomplish the following:
>>> I am solving a coupled system of equations:
>>>
>>> F1(x1,x2,x3) = 0
>>>
>>> F2(x1,x2,x3) = 0
>>>
>>> F3(x1,x2,x3) = 0
>>>
>>>
>>> I am using DMPlex to handle to mesh data. I would like to compute the
>>> norm of residual of each function.
>>>
>>> I assemble the function element wise. I am looking for some interface
>>> like,
>>> ComputNorm(dm, F, fieldIndex).
>>
>>
>> Our thinking is generally that you would use
>>
>>   DMCreateSubDM()
>>
>> and then use VecScatterCreate() with the IS it returns. This will map F into
>> one of the F1/F2/F3.
>>
>> Does this make sense?
>>
>>    Matt
>>
>>>
>>> Thanks
>>> Reddy
>>
>>
>>
>>
>> --
>> 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
>
>
>
> --
> -----------------------------------------------------
> Dharmendar Reddy Palle
> Graduate Student
> Microelectronics Research center,
> University of Texas at Austin,
> 10100 Burnet Road, Bldg. 160
> MER 2.608F, TX 78758-4445
> e-mail: dharmareddy84 at gmail.com
> Phone: +1-512-350-9082
> United States of America.
> Homepage: https://webspace.utexas.edu/~dpr342



-- 
-----------------------------------------------------
Dharmendar Reddy Palle
Graduate Student
Microelectronics Research center,
University of Texas at Austin,
10100 Burnet Road, Bldg. 160
MER 2.608F, TX 78758-4445
e-mail: dharmareddy84 at gmail.com
Phone: +1-512-350-9082
United States of America.
Homepage: https://webspace.utexas.edu/~dpr342