<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"><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"><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" class="gmail_signature" data-smartmail="gmail_signature"><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">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>