<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class=""><br class=""><div><br class=""><blockquote type="cite" class=""><div class="">On Aug 21, 2020, at 8:35 AM, Thibault Bridel-Bertomeu <<a href="mailto:thibault.bridelbertomeu@gmail.com" class="">thibault.bridelbertomeu@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><br clear="all" class=""></div><br class=""><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" class="">knepley@gmail.com</a>> a écrit :<br class=""></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" class=""><div dir="ltr" class="">On Fri, Aug 21, 2020 at 9:10 AM Thibault Bridel-Bertomeu <<a href="mailto:thibault.bridelbertomeu@gmail.com" target="_blank" class="">thibault.bridelbertomeu@gmail.com</a>> wrote:<br class=""></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" class="">Sorry, I sent too soon, I hit the wrong key.<br clear="all" class=""><div class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><br class=""></div><div class="">I wanted to say that context.npoints is the local number of cells.</div><div class=""><br class=""></div><div class="">PetscRHSFunctionImpl allows to generate the hyperbolic part of the right hand side.</div><div class="">Then we have :<br class=""></div><div class=""><br class=""></div><div class=""><div style="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;" class=""><div class=""><span style="color:rgb(32,32,32)" class="">PetscErrorCode</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">PetscIJacobian</span><span style="" class="">(</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">TS</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ts</span><span style="" class="">, </span><span style="color:rgb(160,160,160);font-style:italic" class="">/*!< Time stepping object (see PETSc TS)*/</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">PetscReal</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">t</span><span style="" class="">, </span><span style="color:rgb(160,160,160);font-style:italic" class="">/*!< Current time */</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">Vec</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">Y</span><span style="" class="">, </span><span style="color:rgb(160,160,160);font-style:italic" class="">/*!< Solution vector */</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">Vec</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">Ydot</span><span style="" class="">, </span><span style="color:rgb(160,160,160);font-style:italic" class="">/*!< Time-derivative of solution vector */</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">PetscReal</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">a</span><span style="" class="">, </span><span style="color:rgb(160,160,160);font-style:italic" class="">/*!< Shift */</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">Mat</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">A</span><span style="" class="">, </span><span style="color:rgb(160,160,160);font-style:italic" class="">/*!< Jacobian matrix */</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">Mat</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">B</span><span style="" class="">, </span><span style="color:rgb(160,160,160);font-style:italic" class="">/*!< Preconditioning matrix */</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(9,30,66);font-weight:bold" class="">void</span><span style="" class=""> *</span><span style="color:rgb(32,32,32)" class="">ctxt</span><span style="" class=""> </span><span style="color:rgb(160,160,160);font-style:italic" class="">/*!< Application context */</span></div><div class=""><span style="" class=""> )</span></div><div class=""><span style="" class="">{</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">PETScContext</span><span style="" class=""> *</span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class=""> = (</span><span style="color:rgb(32,32,32)" class="">PETScContext</span><span style="" class="">*) </span><span style="color:rgb(32,32,32)" class="">ctxt</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">HyPar</span><span style="" class=""> *</span><span style="color:rgb(32,32,32)" class="">solver</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">solver</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">_DECLARE_IERR_</span><span style="" class="">;</span></div><br class=""><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">PetscFunctionBegin</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">solver</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">count_IJacobian</span><span style="" class="">++;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">shift</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">a</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">waqt</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">t</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(160,160,160);font-style:italic" class="">/* Construct preconditioning matrix */</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(9,30,66);font-weight:bold" class="">if</span><span style="" class=""> (</span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">flag_use_precon</span><span style="" class="">) { </span><span style="color:rgb(32,32,32)" class="">IERR</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">PetscComputePreconMatImpl</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">B</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">Y</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHECKERR</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">); }</span></div><br class=""><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">PetscFunctionReturn</span><span style="" class="">(</span><span style="color:rgb(101,84,192)" class="">0</span><span style="" class="">);</span></div><div class=""><span style="" class="">}</span></div></div></div></div></div></div></div></div></div><div class=""><br class=""></div><div class="">and<span style="color:rgb(32,32,32)" class=""> 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 class=""></span></div><div class=""><span style="color:rgb(32,32,32)" class="">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 class=""></span></div><div class=""><span style="color:rgb(32,32,32)" class=""><br class=""></span></div><div class=""><div class=""><span class="">/* diagonal element */</span></div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00123" class=""></a><span class=""> </span> <span class="">for</span> (v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pg + v; }</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00124" class=""></a><span class=""> </span> ierr = solver-><a href="http://hypar.github.io/a00144.html#a25c87c7c874c7e08981688fe6c51167b" target="_blank" class="">JFunction</a>(values,(u+nvars*p),solver-><a href="http://hypar.github.io/a00144.html#a48e81a1806b5774943fd9a26e9a190f2" target="_blank" class="">physics</a>,dir,0);</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00125" class=""></a><span class=""> </span> <a href="http://hypar.github.io/a00165.html#a1dab5f429ccd03c4a7fbb605d8ef75b6" target="_blank" class="">_ArrayScale1D_</a>(values,(dxinv*iblank),(nvars*nvars));</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00126" class=""></a><span class=""> </span> ierr = MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES); CHKERRQ(ierr);</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00127" class=""></a><span class=""> </span> <br class=""></div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00128" class=""></a><span class=""> </span> <span class="">/* left neighbor */</span></div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00129" class=""></a><span class=""> </span> <span class="">if</span> (pgL >= 0) {</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00130" class=""></a><span class=""> </span> <span class="">for</span> (v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgL + v; }</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00131" class=""></a><span class=""> </span> ierr = solver-><a href="http://hypar.github.io/a00144.html#a25c87c7c874c7e08981688fe6c51167b" target="_blank" class="">JFunction</a>(values,(u+nvars*pL),solver-><a href="http://hypar.github.io/a00144.html#a48e81a1806b5774943fd9a26e9a190f2" target="_blank" class="">physics</a>,dir,1);</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00132" class=""></a><span class=""> </span> <a href="http://hypar.github.io/a00165.html#a1dab5f429ccd03c4a7fbb605d8ef75b6" target="_blank" class="">_ArrayScale1D_</a>(values,(-dxinv*iblank),(nvars*nvars));</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00133" class=""></a><span class=""> </span> ierr = MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES); CHKERRQ(ierr);</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00134" class=""></a><span class=""> </span> }</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00135" class=""></a><span class=""> </span> </div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00136" class=""></a><span class=""> </span> <span class="">/* right neighbor */</span></div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00137" class=""></a><span class=""> </span> <span class="">if</span> (pgR >= 0) {</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00138" class=""></a><span class=""> </span> <span class="">for</span> (v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgR + v; }</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00139" class=""></a><span class=""> </span> ierr = solver-><a href="http://hypar.github.io/a00144.html#a25c87c7c874c7e08981688fe6c51167b" target="_blank" class="">JFunction</a>(values,(u+nvars*pR),solver-><a href="http://hypar.github.io/a00144.html#a48e81a1806b5774943fd9a26e9a190f2" target="_blank" class="">physics</a>,dir,-1);</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00140" class=""></a><span class=""> </span> <a href="http://hypar.github.io/a00165.html#a1dab5f429ccd03c4a7fbb605d8ef75b6" target="_blank" class="">_ArrayScale1D_</a>(values,(-dxinv*iblank),(nvars*nvars));</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00141" class=""></a><span class=""> </span> ierr = MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES); CHKERRQ(ierr);</div>
<div class=""><a name="m_-549060344189013999_m_-8102338331154726675_l00142" class=""></a><span class=""> </span> }</div><span style="color:rgb(32,32,32)" class=""></span></div><div class=""><span style="color:rgb(32,32,32)" class=""><br class=""></span></div><div class=""><div style="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;" class=""><br class=""></div><div style="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;" class=""><br class=""></div><div style="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;" class="">I do not know if I am clear here ...</div><div style="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;" class="">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 class=""><br class=""></div><div class="">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 class="">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 class=""><br class=""></div><div class="">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></div></div></blockquote><div><br class=""></div><div> Thibault,</div><div><br class=""></div> Since the TS implicit methods end up using SNES internally the option should be available to you without requiring you to be calling the SNES routines directly</div><div><br class=""></div><div> Once you have finalized your approach and if for the implicit case you always work in the snes mf operator mode you can hardwire </div><div><br class=""></div><div> TSGetSNES(ts,&snes);</div><div> SNESSetUseMatrixFree(snes,PETSC_TRUE,PETSC_FALSE);</div><div><br class=""></div><div> in your code so you don't need to always provide the option -snes-mf-operator </div><div><br class=""></div><div> Barry</div><div><br class=""></div><div><br class=""></div><div> <br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div class="gmail_quote"><div class="">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 class=""></div><div class="">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 class=""><pre width="80" class="">TSSetType(ts, TSBEULER);<br class="">DMTSSetIFunctionLocal(dm, <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/DMPlexTSComputeIFunctionFEM.html#DMPlexTSComputeIFunctionFEM" target="_blank" class="">DMPlexTSComputeIFunctionFEM</a>, &ctx);<br class="">DMTSSetIJacobianLocal(dm, <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/DMPlexTSComputeIJacobianFEM.html#DMPlexTSComputeIJacobianFEM" target="_blank" class="">DMPlexTSComputeIJacobianFEM</a>, &ctx);</pre></div><div class=""> with my FV discretization nothing happened, right ?</div><div class=""><br class=""></div><div class="">Thank you, <br class=""></div><div class=""><br class=""></div><div class="">Thibault<br class=""></div><div class=""><br class=""></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" class=""><div class="gmail_quote"><div class=""></div><div class=""> Thanks,</div><div class=""><br class=""></div><div class=""> Matt</div><div class=""> </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" class="">thibault.bridelbertomeu@gmail.com</a>> a écrit :<br class=""></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" class=""><div class="">Hi,</div><div class=""><br class=""></div><div class="">Thanks Matthew and Jed for your input.</div><div class="">I indeed envision an implicit solver in the sense Jed mentioned - Jiri Blazek's book is a nice intro to this concept.</div><div class=""><br class=""></div><div class="">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 class=""></div><div class="">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 class=""><br class=""></div><div class="">Here is the top method for the structured code :</div><div class=""><br class=""></div><div class=""><div style="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;" class=""><div style="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;" class=""><div class=""><span style="color:rgb(9,30,66);font-weight:bold" class="">int</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">total_size</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">.</span><span style="color:rgb(32,32,32)" class="">npoints</span><span style="" class=""> * </span><span style="color:rgb(32,32,32)" class="">solver</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">nvars</span></div></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">TSSetRHSFunction</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ts</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">PETSC_NULL</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">PetscRHSFunctionImpl</span><span style="" class="">,&</span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">SNES</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">snes</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">KSP</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ksp</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">PC</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">pc</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">SNESType</span><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">snestype</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">TSGetSNES</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ts</span><span style="" class="">,&</span><span style="color:rgb(32,32,32)" class="">snes</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">SNESGetType</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">snes</span><span style="" class="">,&</span><span style="color:rgb(32,32,32)" class="">snestype</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><div class=""><span style="" class=""><br class=""></span><span style="color:rgb(160,160,160);font-style:italic" class=""></span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">flag_mat_a</span><span style="" class=""> = </span><span style="color:rgb(101,84,192)" class="">1</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">MatCreateShell</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">MPI_COMM_WORLD</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">total_size</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">total_size</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">PETSC_DETERMINE</span><span style="" class="">,</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">PETSC_DETERMINE</span><span style="" class="">,&</span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">,&</span><span style="color:rgb(32,32,32)" class="">A</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">.</span><span style="color:rgb(32,32,32)" class="">jfnk_eps</span><span style="" class=""> = </span><span style="color:rgb(101,84,192)" class="">1e-7</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">PetscOptionsGetReal</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">NULL</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">NULL</span><span style="" class="">,</span><span style="color:rgb(191,38,0)" class="">"-jfnk_epsilon"</span><span style="" class="">,&</span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">.</span><span style="color:rgb(32,32,32)" class="">jfnk_eps</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">NULL</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">MatShellSetOperation</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">A</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">MATOP_MULT</span><span style="" class="">,(</span><span style="color:rgb(9,30,66);font-weight:bold" class="">void</span><span style="" class=""> (*)(</span><span style="color:rgb(9,30,66);font-weight:bold" class="">void</span><span style="" class="">))</span><span style="color:rgb(32,32,32)" class="">PetscJacobianFunction_JFNK</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">MatSetUp</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">A</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><br class=""><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">.</span><span style="color:rgb(32,32,32)" class="">flag_use_precon</span><span style="" class=""> = </span><span style="color:rgb(101,84,192)" class="">0</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">PetscOptionsGetBool</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">PETSC_NULL</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">PETSC_NULL</span><span style="" class="">,</span><span style="color:rgb(191,38,0)" class="">"-with_pc"</span><span style="" class="">,(</span><span style="color:rgb(32,32,32)" class="">PetscBool</span><span style="" class="">*)(&</span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">.</span><span style="color:rgb(32,32,32)" class="">flag_use_precon</span><span style="" class="">),</span><span style="color:rgb(32,32,32)" class="">PETSC_NULL</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><br class=""><div class=""><span style="" class=""> </span><span style="color:rgb(160,160,160);font-style:italic" class="">/* Set up preconditioner matrix */</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">flag_mat_b</span><span style="" class=""> = </span><span style="color:rgb(101,84,192)" class="">1</span><span style="" class="">;</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">MatCreateAIJ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">MPI_COMM_WORLD</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">total_size</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">total_size</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">PETSC_DETERMINE</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">PETSC_DETERMINE</span><span style="" class="">,</span></div><div class=""><span style="" class=""> (</span><span style="color:rgb(32,32,32)" class="">solver</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">ndims</span><span style="" class="">*</span><span style="color:rgb(101,84,192)" class="">2</span><span style="" class="">+</span><span style="color:rgb(101,84,192)" class="">1</span><span style="" class="">)*</span><span style="color:rgb(32,32,32)" class="">solver</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">nvars</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">NULL</span><span style="" class="">,</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(101,84,192)" class="">2</span><span style="" class="">*</span><span style="color:rgb(32,32,32)" class="">solver</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">ndims</span><span style="" class="">*</span><span style="color:rgb(32,32,32)" class="">solver</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">nvars</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">NULL</span><span style="" class="">,&</span><span style="color:rgb(32,32,32)" class="">B</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">MatSetBlockSize</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">B</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">solver</span><span style="" class="">-></span><span style="color:rgb(32,32,32)" class="">nvars</span><span style="" class="">);</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(160,160,160);font-style:italic" class="">/* Set the RHSJacobian function for TS */</span></div><div class=""><span style="" class=""> </span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class=""> = </span><span style="color:rgb(32,32,32)" class="">TSSetIJacobian</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ts</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">A</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">B</span><span style="" class="">,</span><span style="color:rgb(32,32,32)" class="">PetscIJacobian</span><span style="" class="">,&</span><span style="color:rgb(32,32,32)" class="">context</span><span style="" class="">); </span><span style="color:rgb(32,32,32)" class="">CHKERRQ</span><span style="" class="">(</span><span style="color:rgb(32,32,32)" class="">ierr</span><span style="" class="">);</span></div><div class=""><span style="" class=""><br class=""></span></div></div></div><div class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><div class=""><div class=""><div class=""><div class="">Thibault Bridel-Bertomeu<br class="">—<br class=""></div></div></div></div>Eng, MSc, PhD</div><div class="">Research Engineer</div><div class="">CEA/CESTA</div><div class="">33114 LE BARP</div><div class="">Tel.: (+33)557046924</div><div class="">Mob.: (+33)611025322<br class=""></div><div class="">Mail: <a href="mailto:thibault.bridelbertomeu@gmail.com" target="_blank" class="">thibault.bridelbertomeu@gmail.com</a><br class=""></div></div></div></div></div></div></div><br class=""></div></div><br class=""><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" class="">jed@jedbrown.org</a>> a écrit :<br class=""></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" class="">knepley@gmail.com</a>> writes:<br class="">
<br class="">
> I could never get the FVM stuff to make sense to me for implicit methods.<br class="">
> Here is my problem understanding. If you have an FVM method, it decides<br class="">
> to move "stuff" from one cell to its neighboring cells depending on the<br class="">
> solution to the Riemann problem on each face, which computed the flux. This<br class="">
> is<br class="">
> fine unless the timestep is so big that material can flow through into the<br class="">
> cells beyond the neighbor. Then I should have considered the effect of the<br class="">
> Riemann problem for those interfaces. That would be in the Jacobian, but I<br class="">
> don't know how to compute that Jacobian. I guess you could do everything<br class="">
> matrix-free, but without a preconditioner it seems hard.<br class="">
<br class="">
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 class="">
<br class="">
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 class="">
<br class="">
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 class="">
</blockquote></div>
</blockquote></div>
</blockquote></div><br clear="all" class=""><div class=""><br class=""></div>-- <br class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class="">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">-- Norbert Wiener</div><div class=""><br class=""></div><div class=""><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank" class="">https://www.cse.buffalo.edu/~knepley/</a><br class=""></div></div></div></div></div></div></div></div>
</blockquote></div></div>
</div></blockquote></div><br class=""></body></html>