[petsc-users] matshell for spectral methods in fortran
Luv Sharma
luvsharma11 at gmail.com
Tue May 23 14:36:12 CDT 2017
Dear Barry,
Thanks for your quick reply.
I have tried following:
!--------------------------------------------------------------------------------------------------
module spectral_mech_basic
implicit none
private
#include <petsc/finclude/petsc.h90>
! *PETSc data here*
DM ..
SNES ..
..
contains
!--------------------------------------------------------------------------------------------------
subroutine basicPETSc_init
external :: &
*petsc functions here*
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, &
DMDA_STENCIL_BOX, &
grid(1),grid(2),grid(3), &
1 , 1, worldsize, &
9, 0, &
grid(1),grid(2),localK, &
da,ierr)
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,PETSC_NULL_OBJECT,ierr)
CHKERRQ(ierr)
!call DMCreateMatrix(da,J_shell,ierr)
!CHKERRQ(ierr)
!call DMSetMatType(da,MATSHELL,ierr)
!CHKERRQ(ierr)
!call DMSNESSetJacobianLocal(da,SPEC_mech_formJacobian,PETSC_NULL_OBJECT,ierr) !< function to evaluate stiffness matrix
!CHKERRQ(ierr)
matsize = 9_pInt*grid(1)*grid(2)*grid(3)
call MatCreateShell( PETSC_COMM_WORLD, matsize, matsize, matsize, matsize, PETSC_NULL_OBJECT, J_shell, ierr )
CHKERRQ(ierr)
call SNESSetJacobian( snes, J_shell, J_shell, SPEC_mech_formJacobian, PETSC_NULL_OBJECT, ierr)
CHKERRQ(ierr)
call MatShellSetOperation(J_shell, MATOP_MULT, jac_shell, ierr )
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
call SNESGetKSP(snes,ksp,ierr); CHKERRQ(ierr)
call KSPGetPC(ksp,pc,ierr); CHKERRQ(ierr)
call PCSetType(pc,PCNONE,ierr); CHKERRQ(ierr)
call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr)
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine basicPETSc_init
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function &
basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
implicit none
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
SNESSolve, &
! solve BVP
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr)
CHKERRQ(ierr)
end function BasicPETSc_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the basic residual vector
!--------------------------------------------------------------------------------------------------
subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension(3,3, &
XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &
x_scal
PetscScalar, dimension(3,3, &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal
! constructing residual
…..
….
f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
end subroutine BasicPETSc_formResidual
!—————————————————————————————————————————————————
!> @brief matmult routine
!--------------------------------------------------------------------------------------------------
! a shell jacobian; returns the action J*v
subroutine jac_shell(Jshell,v_in,v_out)
use math, only: &
math_rotate_backward33, &
math_transpose33, &
math_mul3333xx33
use mesh, only: &
grid, &
grid3
use spectral_utilities, only: &
wgt, &
tensorField_real, &
utilities_FFTtensorForward, &
utilities_fourierGammaConvolution, &
utilities_FFTtensorBackward, &
Utilities_constitutiveResponse, &
Utilities_divergenceRMS
implicit none
! DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
! in
PetscScalar, dimension(3,3, &
1000,1,1), intent(in) :: &
v_in
PetscScalar, dimension(3,3, &
1000,1,1), intent(out) :: &
v_out
Mat :: Jshell
PetscErrorCode :: ierr
integer(pInt) :: &
i,j,k,e
e = 0_pInt
tensorField_real = 0.0_pReal
print*, SHAPE(v_in)
print*, SHAPE(v_out)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
e = e + 1_pInt
tensorField_real(1:3,1:3,i,j,k) = j*v
enddo; enddo; enddo
call utilities_FFTtensorForward()
call utilities_fourierGammaConvolution()
call utilities_FFTtensorBackward()
v_out = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
end subroutine jac_shell
subroutine SPEC_mech_formJacobian(snes,xx_local,Jac_pre,Jac,dummy,ierr)
implicit none
SNES :: snes
DM :: dm_local
Vec :: x_local, xx_local
Mat :: Jac_pre, Jac
PetscObject :: dummy
PetscErrorCode :: ierr
end subroutine SPEC_mech_formJacobian
end module spectral_mech_basic
!--------------------------------------------------------------------------------------------------
Best regards,
Luv
> On 23 May 2017, at 21:16, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>
> You didn't include any code related to creating or setting the MATSHELL. What goes wrong with your Fortran code.
>
>
>> On May 23, 2017, at 2:05 PM, Luv Sharma <luvsharma11 at gmail.com> wrote:
>>
>> Dear PETSc team,
>>
>> I am working on a code which solves mechanical equilibrium using spectral methods.
>> I want to make use of the matshell to get the action J*v.
>>
>> I have been able to successfully implement it using petsc4py. But having difficulties to get it working in a fortran code.
>> I am using petsc-3.7.6.
>>
>> Below is a stripped down version of the existing fortran code (module). Can you please help me in figuring out how the right way to do it a code with following structure?
>>
>> !--------------------------------------------------------------------------------------------------
>> module spectral_mech_basic
>>
>> implicit none
>> private
>> #include <petsc/finclude/petsc.h90>
>>
>> ! *PETSc data here*
>> DM ..
>> SNES ..
>> ..
>> contains
>>
>> !--------------------------------------------------------------------------------------------------
>> subroutine basicPETSc_init
>>
>> external :: &
>> *petsc functions here*
>>
>> ! initialize solver specific parts of PETSc
>> call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
>> call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
>> call DMDACreate3d(PETSC_COMM_WORLD, &
>> DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, &
>> DMDA_STENCIL_BOX, &
>> grid(1),grid(2),grid(3), &
>> 1 , 1, worldsize, &
>> 9, 0, &
>> grid(1),grid(2),localK, &
>> da,ierr)
>> CHKERRQ(ierr)
>> call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
>> call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
>> call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,PETSC_NULL_OBJECT,ierr)
>> CHKERRQ(ierr)
>> call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
>> call SNESGetKSP(snes,ksp,ierr); CHKERRQ(ierr)
>> call KSPGetPC(ksp,pc,ierr); CHKERRQ(ierr)
>> call PCSetType(pc,PCNONE,ierr); CHKERRQ(ierr)
>> call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr)
>>
>> ! init fields
>> call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
>> call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
>>
>> end subroutine basicPETSc_init
>> !--------------------------------------------------------------------------------------------------
>>
>> type(tSolutionState) function &
>> basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
>> implicit none
>> ! PETSc Data
>> PetscErrorCode :: ierr
>> SNESConvergedReason :: reason
>> external :: &
>> SNESSolve, &
>> ! solve BVP
>> call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr)
>> CHKERRQ(ierr)
>> end function BasicPETSc_solution
>> !--------------------------------------------------------------------------------------------------
>> !> @brief forms the basic residual vector
>> !--------------------------------------------------------------------------------------------------
>> subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
>>
>> implicit none
>> DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
>> in
>> PetscScalar, dimension(3,3, &
>> XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &
>> x_scal
>> PetscScalar, dimension(3,3, &
>> X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
>> f_scal
>> ! constructing residual
>> …..
>> ….
>>
>> f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
>>
>> end subroutine BasicPETSc_formResidual
>> !--------------------------------------------------------------------------------------------------
>>
>> end module spectral_mech_basic
>> !--------------------------------------------------------------------------------------------------
>>
>> Best regards,
>> Luv
>>
>>> On 3 Nov 2016, at 01:17, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>>
>>>
>>> Is anyone away of cases where PETSc has been used with spectral methods?
>>>
>>> Thanks
>>>
>>> Barry
>>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170523/8ee73acf/attachment-0001.html>
More information about the petsc-users
mailing list