<div dir="ltr"><div dir="ltr"><br clear="all"></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Le ven. 21 août 2020 à 15:23, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> a écrit :<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Fri, Aug 21, 2020 at 9:10 AM Thibault Bridel-Bertomeu <<a href="mailto:thibault.bridelbertomeu@gmail.com" target="_blank">thibault.bridelbertomeu@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Sorry, I sent too soon, I hit the wrong key.<br clear="all"><div><div><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div dir="ltr"><br></div><div>I wanted to say that context.npoints is the local number of cells.</div><div><br></div><div>PetscRHSFunctionImpl allows to generate the hyperbolic part of the right hand side.</div><div>Then we have :<br></div><div><br></div><div><div style="color:rgb(0,0,0);background-color:rgb(255,255,254);font-family:SFMono-Medium,"SF Mono","Segoe UI Mono","Roboto Mono","Ubuntu Mono",Menlo,monospace;font-weight:normal;font-size:13px;line-height:18px;white-space:pre-wrap"><div><span style="color:rgb(32,32,32)">PetscErrorCode</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">PetscIJacobian</span><span style="color:rgb(0,0,0)">(</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">TS</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ts</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(160,160,160);font-style:italic">/*!< Time stepping object (see PETSc TS)*/</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">PetscReal</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">t</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(160,160,160);font-style:italic">/*!< Current time */</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">Vec</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">Y</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(160,160,160);font-style:italic">/*!< Solution vector */</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">Vec</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">Ydot</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(160,160,160);font-style:italic">/*!< Time-derivative of solution vector */</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">PetscReal</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">a</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(160,160,160);font-style:italic">/*!< Shift */</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">Mat</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">A</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(160,160,160);font-style:italic">/*!< Jacobian matrix */</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">Mat</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">B</span><span style="color:rgb(0,0,0)">, </span><span style="color:rgb(160,160,160);font-style:italic">/*!< Preconditioning matrix */</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(9,30,66);font-weight:bold">void</span><span style="color:rgb(0,0,0)"> *</span><span style="color:rgb(32,32,32)">ctxt</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(160,160,160);font-style:italic">/*!< Application context */</span></div><div><span style="color:rgb(0,0,0)"> )</span></div><div><span style="color:rgb(0,0,0)">{</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">PETScContext</span><span style="color:rgb(0,0,0)"> *</span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)"> = (</span><span style="color:rgb(32,32,32)">PETScContext</span><span style="color:rgb(0,0,0)">*) </span><span style="color:rgb(32,32,32)">ctxt</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">HyPar</span><span style="color:rgb(0,0,0)"> *</span><span style="color:rgb(32,32,32)">solver</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">solver</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">_DECLARE_IERR_</span><span style="color:rgb(0,0,0)">;</span></div><br><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">PetscFunctionBegin</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">solver</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">count_IJacobian</span><span style="color:rgb(0,0,0)">++;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">shift</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">a</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">waqt</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">t</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(160,160,160);font-style:italic">/* Construct preconditioning matrix */</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(9,30,66);font-weight:bold">if</span><span style="color:rgb(0,0,0)"> (</span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">flag_use_precon</span><span style="color:rgb(0,0,0)">) { </span><span style="color:rgb(32,32,32)">IERR</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">PetscComputePreconMatImpl</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">B</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">Y</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHECKERR</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">); }</span></div><br><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">PetscFunctionReturn</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(101,84,192)">0</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)">}</span></div></div></div></div></div></div></div></div></div><div><br></div><div>and<span style="color:rgb(32,32,32)"> PetscJacobianFunction_JFNK which I bind to the matrix shell, computes the action of the jacobian on a vector : say U0 is the state of reference and Y the vector upon which to apply the JFNK method, then the PetscJacobianFunction_JFNK returns shift * Y - 1/epsilon * (F(U0 + epsilon*Y) - F(U0)) where F allows to evaluate the hyperbolic flux (shift comes from the TS).<br></span></div><div><span style="color:rgb(32,32,32)">The preconditioning matrix I compute as an approximation to the actual jacobian, that is shift * Identity - Derivative(dF/dU) where dF/dU is, in each cell, a 4x4 matrix that is known exactly for the system of equations I am solving, i.e. Euler equations. For the structured grid, I can loop on the cells and do that 'Derivative' thing at first order by simply taking a finite-difference like approximation with the neighboring cells, Derivative(phi) = phi_i - phi_{i-1} and I assemble the B matrix block by block (JFunction is the dF/dU)<br></span></div><div><span style="color:rgb(32,32,32)"><br></span></div><div><div><span>/* diagonal element */</span></div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00123"></a><span> </span> <span>for</span> (v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pg + v; }</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00124"></a><span> </span> ierr = solver-><a href="http://hypar.github.io/a00144.html#a25c87c7c874c7e08981688fe6c51167b" target="_blank">JFunction</a>(values,(u+nvars*p),solver-><a href="http://hypar.github.io/a00144.html#a48e81a1806b5774943fd9a26e9a190f2" target="_blank">physics</a>,dir,0);</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00125"></a><span> </span> <a href="http://hypar.github.io/a00165.html#a1dab5f429ccd03c4a7fbb605d8ef75b6" target="_blank">_ArrayScale1D_</a>(values,(dxinv*iblank),(nvars*nvars));</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00126"></a><span> </span> ierr = MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES); CHKERRQ(ierr);</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00127"></a><span> </span> <br></div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00128"></a><span> </span> <span>/* left neighbor */</span></div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00129"></a><span> </span> <span>if</span> (pgL >= 0) {</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00130"></a><span> </span> <span>for</span> (v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgL + v; }</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00131"></a><span> </span> ierr = solver-><a href="http://hypar.github.io/a00144.html#a25c87c7c874c7e08981688fe6c51167b" target="_blank">JFunction</a>(values,(u+nvars*pL),solver-><a href="http://hypar.github.io/a00144.html#a48e81a1806b5774943fd9a26e9a190f2" target="_blank">physics</a>,dir,1);</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00132"></a><span> </span> <a href="http://hypar.github.io/a00165.html#a1dab5f429ccd03c4a7fbb605d8ef75b6" target="_blank">_ArrayScale1D_</a>(values,(-dxinv*iblank),(nvars*nvars));</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00133"></a><span> </span> ierr = MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES); CHKERRQ(ierr);</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00134"></a><span> </span> }</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00135"></a><span> </span> </div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00136"></a><span> </span> <span>/* right neighbor */</span></div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00137"></a><span> </span> <span>if</span> (pgR >= 0) {</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00138"></a><span> </span> <span>for</span> (v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgR + v; }</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00139"></a><span> </span> ierr = solver-><a href="http://hypar.github.io/a00144.html#a25c87c7c874c7e08981688fe6c51167b" target="_blank">JFunction</a>(values,(u+nvars*pR),solver-><a href="http://hypar.github.io/a00144.html#a48e81a1806b5774943fd9a26e9a190f2" target="_blank">physics</a>,dir,-1);</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00140"></a><span> </span> <a href="http://hypar.github.io/a00165.html#a1dab5f429ccd03c4a7fbb605d8ef75b6" target="_blank">_ArrayScale1D_</a>(values,(-dxinv*iblank),(nvars*nvars));</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00141"></a><span> </span> ierr = MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES); CHKERRQ(ierr);</div>
<div><a name="m_-549060344189013999_m_-8102338331154726675_l00142"></a><span> </span> }</div><span style="color:rgb(32,32,32)"></span></div><div><span style="color:rgb(32,32,32)"><br></span></div><div><div style="color:rgb(0,0,0);background-color:rgb(255,255,254);font-family:SFMono-Medium,"SF Mono","Segoe UI Mono","Roboto Mono","Ubuntu Mono",Menlo,monospace;font-weight:normal;font-size:13px;line-height:18px;white-space:pre-wrap"><br></div><div style="color:rgb(0,0,0);background-color:rgb(255,255,254);font-family:SFMono-Medium,"SF Mono","Segoe UI Mono","Roboto Mono","Ubuntu Mono",Menlo,monospace;font-weight:normal;font-size:13px;line-height:18px;white-space:pre-wrap"><br></div><div style="color:rgb(0,0,0);background-color:rgb(255,255,254);font-family:SFMono-Medium,"SF Mono","Segoe UI Mono","Roboto Mono","Ubuntu Mono",Menlo,monospace;font-weight:normal;font-size:13px;line-height:18px;white-space:pre-wrap">I do not know if I am clear here ...</div><div style="color:rgb(0,0,0);background-color:rgb(255,255,254);font-family:SFMono-Medium,"SF Mono","Segoe UI Mono","Roboto Mono","Ubuntu Mono",Menlo,monospace;font-weight:normal;font-size:13px;line-height:18px;white-space:pre-wrap">Anyways, I am trying to figure out how to do this shell matrix and this preconditioner using all the FV and DMPlex artillery.</div></div></div></blockquote><div><br></div><div>Okay, that is very clear. We should be able to get the JFNK just with -snes_mf_operator, and put the approximate J construction in DMPlexComputeJacobian_Internal().</div><div>There is an FV section already, and we could just add this. I would need to understand those entries in the pointwise Riemann sense that the other stuff is now.</div></div></div></blockquote><div><br></div><div>Ok i had a quick look and if I understood correctly it would do the job. Setting the snes-mf-operator flag would mean however that we have to go through SNESSetJacobian to set the jacobian and the preconditioning matrix wouldn't it ?</div><div>There might be calls to the Riemann solver to evaluate the dRHS / dU part yes but maybe it's possible to re-use what was computed for the RHS^n ?<br></div><div>In the FV section the jacobian is set to identity which I missed before, but it could explain why when I used the following :<br><pre width="80">TSSetType(ts, TSBEULER);<br>DMTSSetIFunctionLocal(dm, <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/DMPlexTSComputeIFunctionFEM.html#DMPlexTSComputeIFunctionFEM" target="_blank">DMPlexTSComputeIFunctionFEM</a>, &ctx);<br>DMTSSetIJacobianLocal(dm, <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/DMPlexTSComputeIJacobianFEM.html#DMPlexTSComputeIJacobianFEM" target="_blank">DMPlexTSComputeIJacobianFEM</a>, &ctx);</pre></div><div> with my FV discretization nothing happened, right ?</div><div><br></div><div>Thank you, <br></div><div><br></div><div>Thibault<br></div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Le ven. 21 août 2020 à 14:55, Thibault Bridel-Bertomeu <<a href="mailto:thibault.bridelbertomeu@gmail.com" target="_blank">thibault.bridelbertomeu@gmail.com</a>> a écrit :<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi,</div><div><br></div><div>Thanks Matthew and Jed for your input.</div><div>I indeed envision an implicit solver in the sense Jed mentioned - Jiri Blazek's book is a nice intro to this concept.</div><div><br></div><div>Matthew, I do not know exactly what to change right now because although I understand globally what the DMPlexComputeXXXX_Internal methods do, I cannot say for sure line by line what is happening.<br></div><div>In a structured code, I have a an implicit FVM solver with PETSc but I do not use any of the FV structure, not even a DM - I just use C arrays that I transform to PETSc Vec and Mat and build my IJacobian and my preconditioner and gives all that to a TS and it runs. I cannot figure out how to do it with the FV and the DM and all the underlying "shortcuts" that I want to use.</div><div><br></div><div>Here is the top method for the structured code :</div><div><br></div><div><div style="color:rgb(0,0,0);background-color:rgb(255,255,254);font-family:SFMono-Medium,"SF Mono","Segoe UI Mono","Roboto Mono","Ubuntu Mono",Menlo,monospace;font-weight:normal;font-size:13px;line-height:18px;white-space:pre-wrap"><div style="color:rgb(0,0,0);background-color:rgb(255,255,254);font-family:SFMono-Medium,"SF Mono","Segoe UI Mono","Roboto Mono","Ubuntu Mono",Menlo,monospace;font-weight:normal;font-size:13px;line-height:18px;white-space:pre-wrap"><div><span style="color:rgb(9,30,66);font-weight:bold">int</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">total_size</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">.</span><span style="color:rgb(32,32,32)">npoints</span><span style="color:rgb(0,0,0)"> * </span><span style="color:rgb(32,32,32)">solver</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">nvars</span></div></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">TSSetRHSFunction</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ts</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">PETSC_NULL</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">PetscRHSFunctionImpl</span><span style="color:rgb(0,0,0)">,&</span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">SNES</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">snes</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">KSP</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ksp</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">PC</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">pc</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">SNESType</span><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">snestype</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">TSGetSNES</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ts</span><span style="color:rgb(0,0,0)">,&</span><span style="color:rgb(32,32,32)">snes</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">SNESGetType</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">snes</span><span style="color:rgb(0,0,0)">,&</span><span style="color:rgb(32,32,32)">snestype</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)"><br></span><span style="color:rgb(160,160,160);font-style:italic"></span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">flag_mat_a</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(101,84,192)">1</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">MatCreateShell</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">MPI_COMM_WORLD</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">total_size</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">total_size</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">PETSC_DETERMINE</span><span style="color:rgb(0,0,0)">,</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">PETSC_DETERMINE</span><span style="color:rgb(0,0,0)">,&</span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">,&</span><span style="color:rgb(32,32,32)">A</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">.</span><span style="color:rgb(32,32,32)">jfnk_eps</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(101,84,192)">1e-7</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">PetscOptionsGetReal</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">NULL</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">NULL</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(191,38,0)">"-jfnk_epsilon"</span><span style="color:rgb(0,0,0)">,&</span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">.</span><span style="color:rgb(32,32,32)">jfnk_eps</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">NULL</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">MatShellSetOperation</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">A</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">MATOP_MULT</span><span style="color:rgb(0,0,0)">,(</span><span style="color:rgb(9,30,66);font-weight:bold">void</span><span style="color:rgb(0,0,0)"> (*)(</span><span style="color:rgb(9,30,66);font-weight:bold">void</span><span style="color:rgb(0,0,0)">))</span><span style="color:rgb(32,32,32)">PetscJacobianFunction_JFNK</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">MatSetUp</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">A</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><br><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">.</span><span style="color:rgb(32,32,32)">flag_use_precon</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(101,84,192)">0</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">PetscOptionsGetBool</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">PETSC_NULL</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">PETSC_NULL</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(191,38,0)">"-with_pc"</span><span style="color:rgb(0,0,0)">,(</span><span style="color:rgb(32,32,32)">PetscBool</span><span style="color:rgb(0,0,0)">*)(&</span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">.</span><span style="color:rgb(32,32,32)">flag_use_precon</span><span style="color:rgb(0,0,0)">),</span><span style="color:rgb(32,32,32)">PETSC_NULL</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><br><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(160,160,160);font-style:italic">/* Set up preconditioner matrix */</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">flag_mat_b</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(101,84,192)">1</span><span style="color:rgb(0,0,0)">;</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">MatCreateAIJ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">MPI_COMM_WORLD</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">total_size</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">total_size</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">PETSC_DETERMINE</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">PETSC_DETERMINE</span><span style="color:rgb(0,0,0)">,</span></div><div><span style="color:rgb(0,0,0)"> (</span><span style="color:rgb(32,32,32)">solver</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">ndims</span><span style="color:rgb(0,0,0)">*</span><span style="color:rgb(101,84,192)">2</span><span style="color:rgb(0,0,0)">+</span><span style="color:rgb(101,84,192)">1</span><span style="color:rgb(0,0,0)">)*</span><span style="color:rgb(32,32,32)">solver</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">nvars</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">NULL</span><span style="color:rgb(0,0,0)">,</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(101,84,192)">2</span><span style="color:rgb(0,0,0)">*</span><span style="color:rgb(32,32,32)">solver</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">ndims</span><span style="color:rgb(0,0,0)">*</span><span style="color:rgb(32,32,32)">solver</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">nvars</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">NULL</span><span style="color:rgb(0,0,0)">,&</span><span style="color:rgb(32,32,32)">B</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">MatSetBlockSize</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">B</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">solver</span><span style="color:rgb(0,0,0)">-></span><span style="color:rgb(32,32,32)">nvars</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(160,160,160);font-style:italic">/* Set the RHSJacobian function for TS */</span></div><div><span style="color:rgb(0,0,0)"> </span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)"> = </span><span style="color:rgb(32,32,32)">TSSetIJacobian</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ts</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">A</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">B</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(32,32,32)">PetscIJacobian</span><span style="color:rgb(0,0,0)">,&</span><span style="color:rgb(32,32,32)">context</span><span style="color:rgb(0,0,0)">); </span><span style="color:rgb(32,32,32)">CHKERRQ</span><span style="color:rgb(0,0,0)">(</span><span style="color:rgb(32,32,32)">ierr</span><span style="color:rgb(0,0,0)">);</span></div><div><span style="color:rgb(0,0,0)"><br></span></div></div></div><div><div><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div dir="ltr"><div><div><div><div><div>Thibault Bridel-Bertomeu<br>—<br></div></div></div></div>Eng, MSc, PhD</div><div>Research Engineer</div><div>CEA/CESTA</div><div>33114 LE BARP</div><div>Tel.: (+33)557046924</div><div>Mob.: (+33)611025322<br></div><div>Mail: <a href="mailto:thibault.bridelbertomeu@gmail.com" target="_blank">thibault.bridelbertomeu@gmail.com</a><br></div></div></div></div></div></div></div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Le jeu. 20 août 2020 à 18:43, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> a écrit :<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> writes:<br>
<br>
> I could never get the FVM stuff to make sense to me for implicit methods.<br>
> Here is my problem understanding. If you have an FVM method, it decides<br>
> to move "stuff" from one cell to its neighboring cells depending on the<br>
> solution to the Riemann problem on each face, which computed the flux. This<br>
> is<br>
> fine unless the timestep is so big that material can flow through into the<br>
> cells beyond the neighbor. Then I should have considered the effect of the<br>
> Riemann problem for those interfaces. That would be in the Jacobian, but I<br>
> don't know how to compute that Jacobian. I guess you could do everything<br>
> matrix-free, but without a preconditioner it seems hard.<br>
<br>
So long as we're using method of lines, the flux is just instantaneous flux, not integrated over some time step. It has the same meaning for implicit and explicit.<br>
<br>
An explicit method would be unstable if you took such a large time step (CFL) and an implicit method will not simultaneously be SSP and higher than first order, but it's still a consistent discretization of the problem.<br>
<br>
It's common (done in FUN3D and others) to precondition with a first-order method, where gradient reconstruction/limiting is skipped. That's what I'd recommend because limiting creates nasty nonlinearities and the resulting discretizations lack h-ellipticity which makes them very hard to solve.<br>
</blockquote></div>
</blockquote></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote></div></div>