<html><head><meta http-equiv="Content-Type" content="text/html charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class="">Dear Barry, <div class=""><br class=""></div><div class="">Thanks for your quick reply.</div><div class="">I have tried following:</div><div class=""><br class=""></div><div class=""><div class=""><div class="">!--------------------------------------------------------------------------------------------------</div></div><div class=""><div class="">module spectral_mech_basic</div><div class=""><br class=""></div><div class=""> implicit none</div><div class=""> private</div><div class="">#include <petsc/finclude/petsc.h90></div><div class=""><br class=""></div><div class="">! *PETSc data here*</div><div class=""> DM ..</div><div class=""> SNES ..</div><div class=""> ..</div><div class="">contains</div><div class=""><br class=""></div><div class="">!--------------------------------------------------------------------------------------------------</div><div class="">subroutine basicPETSc_init</div><div class=""><br class=""></div><div class=""> external :: &</div><div class=""> *petsc functions here*</div><div class=""><br class=""></div><div class="">! initialize solver specific parts of PETSc</div><div class=""> call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)</div><div class=""> call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) </div><div class=""> call DMDACreate3d(PETSC_COMM_WORLD, &</div><div class=""> DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & </div><div class=""> DMDA_STENCIL_BOX, & </div><div class=""> grid(1),grid(2),grid(3), & </div><div class=""> 1 , 1, worldsize, &</div><div class=""> 9, 0, & </div><div class=""> grid(1),grid(2),localK, & </div><div class=""> da,ierr) </div><div class=""> CHKERRQ(ierr)</div><div class=""> call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)</div><div class=""> call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) </div><div class=""> call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,PETSC_NULL_OBJECT,ierr) </div><div class=""> CHKERRQ(ierr) </div><div class=""><br class=""></div><div class=""><div class=""><br class=""></div><div class=""><font color="#0433ff" class=""> !call DMCreateMatrix(da,J_shell,ierr)</font></div><div class=""><font color="#0433ff" class=""> !CHKERRQ(ierr)</font></div><div class=""><font color="#0433ff" class=""> !call DMSetMatType(da,MATSHELL,ierr)</font></div><div class=""><font color="#0433ff" class=""> !CHKERRQ(ierr)</font></div><div class=""><div class=""><font color="#0433ff" class=""> !call DMSNESSetJacobianLocal(da,SPEC_mech_formJacobian,PETSC_NULL_OBJECT,ierr) !< function to evaluate stiffness matrix</font></div><div class=""><font color="#0433ff" class=""> !CHKERRQ(ierr)</font></div><div class=""><font color="#0433ff" class=""><br class=""></font></div></div><div class=""><font color="#0433ff" class=""><br class=""></font></div><div class=""><font color="#0433ff" class=""> matsize = 9_pInt*grid(1)*grid(2)*grid(3)</font></div><div class=""><font color="#0433ff" class=""> call MatCreateShell( PETSC_COMM_WORLD, matsize, matsize, matsize, matsize, PETSC_NULL_OBJECT, J_shell, ierr )</font></div><div class=""><font color="#0433ff" class=""> CHKERRQ(ierr)</font></div><div class=""><br class=""></div><div class=""><font color="#0433ff" class=""> call SNESSetJacobian( snes, J_shell, J_shell, SPEC_mech_formJacobian, PETSC_NULL_OBJECT, ierr)</font></div><div class=""><font color="#0433ff" class=""> CHKERRQ(ierr)</font></div><div class=""><font color="#0433ff" class=""><br class=""></font></div><div class=""><font color="#0433ff" class=""> call MatShellSetOperation(J_shell, MATOP_MULT, jac_shell, ierr )</font></div><div class=""><font color="#0433ff" class=""> CHKERRQ(ierr)</font></div></div><div class=""><br class=""></div><div class=""> call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) </div><div class=""> call SNESGetKSP(snes,ksp,ierr); CHKERRQ(ierr) </div><div class=""> call KSPGetPC(ksp,pc,ierr); CHKERRQ(ierr) </div><div class=""> call PCSetType(pc,PCNONE,ierr); CHKERRQ(ierr) </div><div class=""> call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) </div><div class=""><br class=""></div><div class="">! init fields </div><div class=""> call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) </div><div class=""> call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) </div><div class=""><br class=""></div><div class="">end subroutine basicPETSc_init</div><div class="">!--------------------------------------------------------------------------------------------------</div><div class=""><br class=""></div><div class="">type(tSolutionState) function &</div><div class=""> basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)</div><div class=""> implicit none</div><div class="">! PETSc Data</div><div class=""> PetscErrorCode :: ierr </div><div class=""> SNESConvergedReason :: reason</div><div class=""> external :: &</div><div class=""> SNESSolve, &</div><div class="">! solve BVP </div><div class=""><span class="" style="font-size: 14px;"> call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr)</span></div><div class=""> CHKERRQ(ierr)</div><div class="">end function BasicPETSc_solution</div><div class="">!--------------------------------------------------------------------------------------------------</div><div class="">!> @brief forms the basic residual vector</div><div class="">!--------------------------------------------------------------------------------------------------</div><div class="">subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)</div><div class=""><br class=""></div><div class=""> implicit none</div><div class=""> DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &</div><div class=""> in</div><div class=""> PetscScalar, dimension(3,3, &</div><div class=""> XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &</div><div class=""> x_scal</div><div class=""> PetscScalar, dimension(3,3, &</div><div class=""> X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &</div><div class=""> f_scal</div><div class="">! constructing residual</div><div class="">…..</div><div class="">….</div><div class=""><br class=""></div><div class=""> f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)</div><div class=""><br class=""></div><div class="">end subroutine BasicPETSc_formResidual</div><div class="">!—————————————————————————————————————————————————</div><div class=""><div class=""><font color="#0433ff" class="">!> @brief matmult routine</font></div><div class=""><font color="#0433ff" class="">!--------------------------------------------------------------------------------------------------</font></div><div class=""><font color="#0433ff" class=""> ! a shell jacobian; returns the action J*v</font></div><div class=""><font color="#0433ff" class="">subroutine jac_shell(Jshell,v_in,v_out)</font></div><div class=""><font color="#0433ff" class=""><br class=""></font></div><div class=""><font color="#0433ff" class=""> use math, only: &</font></div><div class=""><font color="#0433ff" class=""> math_rotate_backward33, &</font></div><div class=""><font color="#0433ff" class=""> math_transpose33, &</font></div><div class=""><font color="#0433ff" class=""> math_mul3333xx33</font></div><div class=""><font color="#0433ff" class=""> use mesh, only: &</font></div><div class=""><font color="#0433ff" class=""> grid, &</font></div><div class=""><font color="#0433ff" class=""> grid3</font></div><div class=""><font color="#0433ff" class=""> use spectral_utilities, only: &</font></div><div class=""><font color="#0433ff" class=""> wgt, &</font></div><div class=""><font color="#0433ff" class=""> tensorField_real, &</font></div><div class=""><font color="#0433ff" class=""> utilities_FFTtensorForward, &</font></div><div class=""><font color="#0433ff" class=""> utilities_fourierGammaConvolution, &</font></div><div class=""><font color="#0433ff" class=""> utilities_FFTtensorBackward, &</font></div><div class=""><font color="#0433ff" class=""> Utilities_constitutiveResponse, &</font></div><div class=""><font color="#0433ff" class=""> Utilities_divergenceRMS</font></div><div class=""><font color="#0433ff" class=""> implicit none</font></div><div class=""><font color="#0433ff" class="">! DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &</font></div><div class=""><font color="#0433ff" class="">! in</font></div><div class=""><font color="#0433ff" class=""> PetscScalar, dimension(3,3, &</font></div><div class=""><font color="#0433ff" class=""> 1000,1,1), intent(in) :: &</font></div><div class=""><font color="#0433ff" class=""> v_in</font></div><div class=""><font color="#0433ff" class=""> PetscScalar, dimension(3,3, &</font></div><div class=""><font color="#0433ff" class=""> 1000,1,1), intent(out) :: &</font></div><div class=""><font color="#0433ff" class=""> v_out</font></div><div class=""><font color="#0433ff" class=""> </font></div><div class=""><font color="#0433ff" class=""> Mat :: Jshell</font></div><div class=""><font color="#0433ff" class=""> PetscErrorCode :: ierr</font></div><div class=""><font color="#0433ff" class=""> integer(pInt) :: &</font></div><div class=""><font color="#0433ff" class=""> i,j,k,e</font></div><div class=""><font color="#0433ff" class=""><br class=""></font></div><div class=""><font color="#0433ff" class=""> e = 0_pInt</font></div><div class=""><font color="#0433ff" class=""><br class=""></font></div><div class=""><font color="#0433ff" class=""> tensorField_real = 0.0_pReal</font></div><div class=""><font color="#0433ff" class=""> print*, SHAPE(v_in)</font></div><div class=""><font color="#0433ff" class=""> print*, SHAPE(v_out)</font></div><div class=""><font color="#0433ff" class=""> </font></div><div class=""><font color="#0433ff" class=""> do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)</font></div><div class=""><font color="#0433ff" class=""> e = e + 1_pInt</font></div><div class=""><font color="#0433ff" class=""> tensorField_real(1:3,1:3,i,j,k) = j*v</font></div><div class=""><font color="#0433ff" class=""> enddo; enddo; enddo</font></div><div class=""><font color="#0433ff" class=""> call utilities_FFTtensorForward()</font></div><div class=""><font color="#0433ff" class=""> call utilities_fourierGammaConvolution()</font></div><div class=""><font color="#0433ff" class=""> call utilities_FFTtensorBackward()</font></div><div class=""><font color="#0433ff" class=""> v_out = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)</font></div><div class=""><font color="#0433ff" class=""><br class=""></font></div><div class=""><font color="#0433ff" class="">end subroutine jac_shell</font></div><div class=""><font color="#0433ff" class=""><br class=""></font></div><div class=""><font color="#0433ff" class=""><br class=""></font></div><div class=""><font color="#0433ff" class="">subroutine SPEC_mech_formJacobian(snes,xx_local,Jac_pre,Jac,dummy,ierr)</font></div><div class=""><font color="#0433ff" class=""> implicit none</font></div><div class=""><font color="#0433ff" class=""> SNES :: snes</font></div><div class=""><font color="#0433ff" class=""> DM :: dm_local</font></div><div class=""><font color="#0433ff" class=""> Vec :: x_local, xx_local</font></div><div class=""><font color="#0433ff" class=""> Mat :: Jac_pre, Jac</font></div><div class=""><font color="#0433ff" class=""> PetscObject :: dummy</font></div><div class=""><font color="#0433ff" class=""> PetscErrorCode :: ierr</font></div><div class=""><font color="#0433ff" class=""><br class=""></font></div><div class=""><font color="#0433ff" class="">end subroutine SPEC_mech_formJacobian</font></div><div class=""><br class=""></div></div><div class=""><br class=""></div><div class=""><br class=""></div><div class="">end module spectral_mech_basic</div></div><div class="">!--------------------------------------------------------------------------------------------------</div></div><div class=""><br class=""></div><div class="">Best regards,</div><div class="">Luv</div><div class=""><br class=""></div><div class=""><br class=""></div><div class=""><br class=""></div><div class=""><br class=""><div><blockquote type="cite" class=""><div class="">On 23 May 2017, at 21:16, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" class="">bsmith@mcs.anl.gov</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div class=""><br class=""> You didn't include any code related to creating or setting the MATSHELL. What goes wrong with your Fortran code.<br class=""><br class=""><br class=""><blockquote type="cite" class="">On May 23, 2017, at 2:05 PM, Luv Sharma <<a href="mailto:luvsharma11@gmail.com" class="">luvsharma11@gmail.com</a>> wrote:<br class=""><br class="">Dear PETSc team,<br class=""><br class="">I am working on a code which solves mechanical equilibrium using spectral methods.<br class="">I want to make use of the matshell to get the action J*v.<br class=""><br class="">I have been able to successfully implement it using petsc4py. But having difficulties to get it working in a fortran code.<br class="">I am using petsc-3.7.6. <br class=""><br class="">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?<br class=""><br class="">!--------------------------------------------------------------------------------------------------<br class="">module spectral_mech_basic<br class=""><br class=""> implicit none<br class=""> private<br class="">#include <petsc/finclude/petsc.h90><br class=""><br class="">! *PETSc data here*<br class=""> DM ..<br class=""> SNES ..<br class=""> ..<br class="">contains<br class=""><br class="">!--------------------------------------------------------------------------------------------------<br class="">subroutine basicPETSc_init<br class=""><br class=""> external :: &<br class=""> *petsc functions here*<br class=""><br class="">! initialize solver specific parts of PETSc<br class=""> call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)<br class=""> call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) <br class=""> call DMDACreate3d(PETSC_COMM_WORLD, &<br class=""> DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & <br class=""> DMDA_STENCIL_BOX, & <br class=""> grid(1),grid(2),grid(3), & <br class=""> 1 , 1, worldsize, &<br class=""> 9, 0, & <br class=""> grid(1),grid(2),localK, & <br class=""> da,ierr) <br class=""> CHKERRQ(ierr)<br class=""> call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)<br class=""> call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) <br class=""> call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,PETSC_NULL_OBJECT,ierr) <br class=""> CHKERRQ(ierr) <br class=""> call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) <br class=""> call SNESGetKSP(snes,ksp,ierr); CHKERRQ(ierr) <br class=""> call KSPGetPC(ksp,pc,ierr); CHKERRQ(ierr) <br class=""> call PCSetType(pc,PCNONE,ierr); CHKERRQ(ierr) <br class=""> call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) <br class=""><br class="">! init fields <br class=""> call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) <br class=""> call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) <br class=""><br class="">end subroutine basicPETSc_init<br class="">!--------------------------------------------------------------------------------------------------<br class=""><br class="">type(tSolutionState) function &<br class=""> basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)<br class=""> implicit none<br class="">! PETSc Data<br class=""> PetscErrorCode :: ierr <br class=""> SNESConvergedReason :: reason<br class=""> external :: &<br class=""> SNESSolve, &<br class="">! solve BVP <br class=""> call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr)<br class=""> CHKERRQ(ierr)<br class="">end function BasicPETSc_solution<br class="">!--------------------------------------------------------------------------------------------------<br class="">!> @brief forms the basic residual vector<br class="">!--------------------------------------------------------------------------------------------------<br class="">subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)<br class=""><br class=""> implicit none<br class=""> DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &<br class=""> in<br class=""> PetscScalar, dimension(3,3, &<br class=""> XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &<br class=""> x_scal<br class=""> PetscScalar, dimension(3,3, &<br class=""> X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &<br class=""> f_scal<br class="">! constructing residual<br class="">…..<br class="">….<br class=""><br class=""> f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)<br class=""><br class="">end subroutine BasicPETSc_formResidual<br class="">!--------------------------------------------------------------------------------------------------<br class=""><br class="">end module spectral_mech_basic<br class="">!--------------------------------------------------------------------------------------------------<br class=""><br class="">Best regards,<br class="">Luv<br class=""><br class=""><blockquote type="cite" class="">On 3 Nov 2016, at 01:17, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" class="">bsmith@mcs.anl.gov</a>> wrote:<br class=""><br class=""><br class=""> Is anyone away of cases where PETSc has been used with spectral methods?<br class=""><br class=""> Thanks<br class=""><br class=""> Barry<br class=""><br class=""></blockquote></blockquote><br class=""></div></div></blockquote></div><br class=""></div></body></html>