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 <[email protected]> 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 <[email protected]> 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 <[email protected]> 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 <[email protected]> wrote: >>>> >>>> >>>> Is anyone away of cases where PETSc has been used with spectral methods? >>>> >>>> Thanks >>>> >>>> Barry >>>> >> >
