<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=us-ascii"><meta name=Generator content="Microsoft Word 14 (filtered medium)"><style><!--
/* Font Definitions */
@font-face
{font-family:SimSun;
panose-1:2 1 6 0 3 1 1 1 1 1;}
@font-face
{font-family:SimSun;
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:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
@font-face
{font-family:SimSun;
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:12.0pt;
font-family:"Times New Roman","serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.MsoListParagraph, li.MsoListParagraph, div.MsoListParagraph
{mso-style-priority:34;
margin-top:0in;
margin-right:0in;
margin-bottom:0in;
margin-left:.5in;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
span.EmailStyle17
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.25in 1.0in 1.25in;}
div.WordSection1
{page:WordSection1;}
/* List Definitions */
@list l0
{mso-list-id:2047951905;
mso-list-type:hybrid;
mso-list-template-ids:2062306650 67698705 67698713 67698715 67698703 67698713 67698715 67698703 67698713 67698715;}
@list l0:level1
{mso-level-text:"%1\)";
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level2
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level3
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level4
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level5
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level6
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level7
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level8
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level9
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
ol
{margin-bottom:0in;}
ul
{margin-bottom:0in;}
--></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-US link=blue vlink=purple><div class=WordSection1><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'>From:</span></b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'> Matthew Knepley [mailto:knepley@gmail.com] <br><b>Sent:</b> Tuesday, June 04, 2013 12:32 PM<br><b>To:</b> Jin, Shuangshuang<br><b>Cc:</b> petsc-users@mcs.anl.gov<br><b>Subject:</b> Re: [petsc-users] a question about DAE Function<o:p></o:p></span></p><p class=MsoNormal><o:p> </o:p></p><div><p class=MsoNormal>On Tue, Jun 4, 2013 at 3:13 PM, Jin, Shuangshuang <<a href="mailto:Shuangshuang.Jin@pnnl.gov" target="_blank">Shuangshuang.Jin@pnnl.gov</a>> wrote:<o:p></o:p></p><div><div><blockquote style='border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in'><div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>Hello, I have a question regarding how to set up the IFunction for TSSetIFunction. Please see the detailed illustration below:<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'> <o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>Inside my IFunction, I have one of the f equations defined as:<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'> <o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>if (i >= xstart && i < xstart+xlen) {<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'> f[i+2] = sin(x[i])*x[i+2] + cos(x[i])*x[i+3];<o:p></o:p></span></p></div></div></blockquote><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoListParagraph style='text-indent:-.25in;mso-list:l0 level1 lfo1'><![if !supportLists]><span style='mso-list:Ignore'>1)<span style='font:7.0pt "Times New Roman"'> </span></span><![endif]>You should not be setting f[i+2]. You always set the residual local to your process.<o:p></o:p></p><p class=MsoNormal><b><i><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>[Jin, Shuangshuang] </span></i></b><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'> Sorry, </span><span style='color:#1F497D'>I typed it wrong, the condition is “if (i+2 >= xstart && i+2 < xstart+xlen) {“ instead.</span><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p></o:p></span></p></div><div><p class=MsoNormal> <o:p></o:p></p></div><blockquote style='border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in'><div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'> for (k = 0; k < n; k++) {<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'> f[i+2] += -(cos(x[k])*a[k])*b) - (sin(x[k])*c[k])*d); (1)<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'> }<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>}<o:p></o:p></span></p></div></div></blockquote><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>2) This in effect says that your Jacobian is completely dense, or that your problem is fully coupled. I see two<o:p></o:p></p></div><div><p class=MsoNormal> choices:<o:p></o:p></p></div><div><p class=MsoNormal><b><i><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>[Jin, Shuangshuang] I have four sets of f functions: f[i], f[i+1], f[i+2], and f[i+3]. (i=1..n), the Jacobian matrix for the whole DAE problem is sparse, but the problem is fully coupled. <o:p></o:p></span></i></b></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p></div><div><p class=MsoNormal> a) Use scalable dense linear algebra, like Elemental, and formulate your operations as linear algebra<span style='color:#1F497D'><o:p></o:p></span></p></div><div><p class=MsoNormal><b><i><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>[Jin, Shuangshuang] Do you have a sample source code on the web for this?<o:p></o:p></span></i></b></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p></div><div><p class=MsoNormal> b) If this is a boundary integral equation, consider using the FMM technology<o:p></o:p></p><p class=MsoNormal><b><i><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>[Jin, Shuangshuang] It’s not.</span></i></b><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p></o:p></span></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>Either way, this is not a PDE, and so using that PETSc parallel paradigm is wrong.<o:p></o:p></p><p class=MsoNormal><b><i><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>[Jin, Shuangshuang] This is a DAE, I set up the section 6.1.2 Solving Differential Algebraic Equations of the Manual, form IFUNCITON, IJacobain, and set up initial values for x. If I scatter x every time inside IFUNCTION, the problem is solved in parallel when I use multiple threads. I just felt that my way of scattering x is inefficient. I don’t understand why using the PETSc parallel paradigm to solve the problem is wrong?<o:p></o:p></span></i></b></p><p class=MsoNormal><b><i><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></i></b></p><p class=MsoNormal><b><i><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Thanks,<o:p></o:p></span></i></b></p><p class=MsoNormal><b><i><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Shuangshuang</span></i></b><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p></o:p></span></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal> Matt<o:p></o:p></p></div><div><p class=MsoNormal> <o:p></o:p></p></div><blockquote style='border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in'><div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>Where x is an distributed array, a and c are constant arrays, b and d are constants. n is the length of the global length of x.<o:p></o:p></span></p></div><div><p class=MsoNormal style='text-indent:9.0pt'> <span style='font-family:"Calibri","sans-serif"'><o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>The problem is, since x is distributed across all the processes, each process only owns a piece of items of x, when I do (1) as written above, no matter which processor f[i+2] is on, it has to grab all the x array data from other processors, which means I have to do a scattering of x array in advance to make the x[k] value available. This is quite an inefficient way given the scattering of vector would happen in each iteration when the IFuntion is called.<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'> <o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>I am wondering if this is the only way I have to deal with this problem, or am I trapped in a wrong direction?<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'> <o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>Your comments are highly appreciated!<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'> <o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>Thanks,<o:p></o:p></span></p></div><div><p class=MsoNormal><span style='font-family:"Calibri","sans-serif"'>Shuangshuang<o:p></o:p></span></p></div><div><p class=MsoNormal> <span style='font-family:"Calibri","sans-serif"'><o:p></o:p></span></p></div><div><p class=MsoNormal> <span style='font-family:"Calibri","sans-serif"'><o:p></o:p></span></p></div></div></blockquote></div><p class=MsoNormal><br><br clear=all><o:p></o:p></p><div><p class=MsoNormal><o:p> </o:p></p></div><p class=MsoNormal>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener <o:p></o:p></p></div></div></div></body></html>