[petsc-users] matshell for spectral methods in fortran
Barry Smith
bsmith at mcs.anl.gov
Tue May 23 15:20:11 CDT 2017
I cannot easily see why this would or wouldn't work. If you send me the entire code as an attachment and makefile I can try to run it.
Barry
> On May 23, 2017, at 2:36 PM, Luv Sharma <luvsharma11 at gmail.com> wrote:
>
> 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
>>>>
>>
>
More information about the petsc-users
mailing list