<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"><head><meta http-equiv=Content-Type content="text/html; charset=utf-8"><meta name=Generator content="Microsoft Word 15 (filtered medium)"><style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:DengXian;
panose-1:2 1 6 0 3 1 1 1 1 1;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:"\@DengXian";
panose-1:2 1 6 0 3 1 1 1 1 1;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:#0563C1;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:#954F72;
text-decoration:underline;}
pre
{mso-style-priority:99;
mso-style-link:"HTML Preformatted Char";
margin:0in;
margin-bottom:.0001pt;
font-size:10.0pt;
font-family:"Courier New";}
span.EmailStyle17
{mso-style-type:personal-compose;
font-family:"Calibri",sans-serif;
color:windowtext;}
span.HTMLPreformattedChar
{mso-style-name:"HTML Preformatted Char";
mso-style-priority:99;
mso-style-link:"HTML Preformatted";
font-family:"Courier New";}
span.line
{mso-style-name:line;}
span.hljs-keyword
{mso-style-name:hljs-keyword;}
span.hljs-comment
{mso-style-name:hljs-comment;}
span.hljs-number
{mso-style-name:hljs-number;}
span.hljs-builtin
{mso-style-name:hljs-built_in;}
span.hljs-string
{mso-style-name:hljs-string;}
.MsoChpDefault
{mso-style-type:export-only;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]--></head><body lang=en-DE link="#0563C1" vlink="#954F72"><div class=WordSection1><p class=MsoNormal><span lang=EN-US>Dear PETSc Experts,<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US>I am implementing a matrix-free jacobian for my SNES solver in Fortran. (command line option -snes_type newtonls -ksp_type gmres)<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US>In the main program, I define my residual and jacobian and matrix-free jacobian like the following,<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US>…<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US>call DMDASNESSetFunctionLocal(DM_mech, INSERT_VALUES, formResidual, PETSC_NULL_SNES, err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US>call DMSNESSetJacobianLocal(DM_mech, formJacobian, PETSC_NULL_SNES, err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US>…<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'>subroutine formJacobian(residual_subdomain,F,Jac_pre,Jac,dummy,err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> <o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'>#include <petsc/finclude/petscmat.h><o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> use petscmat<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> implicit None<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> residual_subdomain !< DMDA info (needs to be named "in" for macros like XRANGE to work)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: &<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> F !< deformation gradient field<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> Mat :: Jac, Jac_pre<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> PetscObject :: dummy<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> PetscErrorCode :: err_PETSc<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> PetscInt :: N_dof ! global number of DoF, maybe only a placeholder<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> N_dof = 9*product(cells(1:2))*cells3 <o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> print*, 'in my jac'<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> <o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N_dof,N_dof,0,Jac,err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> CHKERRQ(err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> call MatShellSetOperation(Jac,MATOP_MULT,GK_op,err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> CHKERRQ(err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> <o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> print*, 'in my jac'<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> ! for jac preconditioner<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N_dof,N_dof,0,Jac_pre,err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> CHKERRQ(err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> call MatShellSetOperation(Jac_pre,MATOP_MULT,GK_op,err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> CHKERRQ(err_PETSc)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'> print*, 'in my jac'<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:"Courier New"'>end subroutine formJacobian<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p><pre><span class=hljs-keyword>subroutine</span><span class=line> GK_op(Jac,dF,output,err_PETSc)</span><o:p></o:p></pre><pre><span class=line> </span><span class=hljs-keyword>real</span><span class=line>(pREAL), </span><span class=hljs-keyword>dimension</span><span class=line>(</span><span class=hljs-number>3</span><span class=line>,</span><span class=hljs-number>3</span><span class=line>,cells(</span><span class=hljs-number>1</span><span class=line>),cells(</span><span class=hljs-number>2</span><span class=line>),cells3), </span><span class=hljs-keyword>intent</span><span class=line>(</span><span class=hljs-keyword>in</span><span class=line>) :: &</span><o:p></o:p></pre><pre><span class=line> dF </span><o:p></o:p></pre><pre><span class=line> </span><span class=hljs-keyword>real</span><span class=line>(pREAL), </span><span class=hljs-keyword>dimension</span><span class=line>(</span><span class=hljs-number>3</span><span class=line>,</span><span class=hljs-number>3</span><span class=line>,cells(</span><span class=hljs-number>1</span><span class=line>),cells(</span><span class=hljs-number>2</span><span class=line>),cells3), </span><span class=hljs-keyword>intent</span><span class=line>(</span><span class=hljs-keyword>out</span><span class=line>) :: &</span><o:p></o:p></pre><pre><span class=line> output </span><o:p></o:p></pre><pre><span class=line> </span><span class=hljs-keyword>real</span><span class=line>(pREAL), </span><span class=hljs-keyword>dimension</span><span class=line>(</span><span class=hljs-number>3</span><span class=line>,</span><span class=hljs-number>3</span><span class=line>) :: &</span><o:p></o:p></pre><pre><span class=line> deltaF_aim = </span><span class=hljs-number>0.0_pREAL</span><o:p></o:p></pre><pre><o:p> </o:p></pre><pre><span class=line> Mat :: Jac</span><o:p></o:p></pre><pre><span class=line> PetscErrorCode :: err_PETSc</span><o:p></o:p></pre><pre><o:p> </o:p></pre><pre><span class=line> </span><span class=hljs-keyword>integer</span><span class=line> :: i, j, k, e</span><o:p></o:p></pre><pre><o:p> </o:p></pre><pre><span class=line> <span lang=EN-US>… a lot of calculations …<o:p></o:p></span></span></pre><pre><span class=line><span lang=EN-US><o:p> </o:p></span></span></pre><pre><span class=line> </span><span class=hljs-builtin>print</span><span class=line>*, </span><span class=hljs-string>'in GK op'</span><o:p></o:p></pre><pre><span class=line> </span><o:p></o:p></pre><pre><span class=hljs-keyword>end</span><span class=line> </span><span class=hljs-keyword>subroutine</span><span class=line> GK_op</span><o:p></o:p></pre><pre><o:p> </o:p></pre><p class=MsoNormal><span lang=EN-US>The first question is that: it seems I still need to explicitly define the interface of MatCreateShell() and MatShellSetOperation() to properly use them, even though I include them via “use petscmat”. It is a little bit strange to me, since some examples do not perform this step. <o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US>Then the main issue is that I can build my own Jacobian from my call back function formJacobian, and confirm my Jacobian is a shell matrix (by MatView). However, my customized operator GK_op is not called when solving the nonlinear system (not print my “in GK op”). When I try to monitor my SNES, it gives me some conventional output not mentioning my matrix-free operations. So I guess my customized MATOP_MULT may be not associated with Jacobian. Or my configuration is somehow wrong. Could you help me solve this issue? <o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p><p class=MsoNormal><span lang=EN-US>Thanks,<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US>Yi<o:p></o:p></span></p><p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p></div><BR />
<BR />
<HR />
-------------------------------------------------<BR />
Stay up to date and follow us on LinkedIn, Twitter and YouTube.<BR />
<BR />
Max-Planck-Institut für Eisenforschung GmbH<BR />
Max-Planck-Straße 1<BR />
D-40237 Düsseldorf<BR />
<BR />
Handelsregister B 2533 <BR />
Amtsgericht Düsseldorf<BR />
<BR />
Geschäftsführung<BR />
Prof. Dr. Gerhard Dehm<BR />
Prof. Dr. Jörg Neugebauer<BR />
Prof. Dr. Dierk Raabe<BR />
Dr. Kai de Weldige<BR />
<BR />
Ust.-Id.-Nr.: DE 11 93 58 514 <BR />
Steuernummer: 105 5891 1000<BR />
<BR />
<BR />
Please consider that invitations and e-mails of our institute are <BR />
only valid if they end with …@mpie.de. <BR />
If you are not sure of the validity please contact rco@mpie.de<BR />
<BR />
Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails<BR />
aus unserem Haus nur mit der Endung …@mpie.de gültig sind. <BR />
In Zweifelsfällen wenden Sie sich bitte an rco@mpie.de<BR />
-------------------------------------------------<BR />
</body></html>