<div dir="ltr"><div dir="ltr"><div>Thanks Matthew for your very fast answer. Sorry I sent another mail at the same time my first one wasn't complete.</div><div><br></div></div><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Le ven. 21 août 2020 à 15:09, 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 8:56 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"><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></blockquote><div><br></div><div>Unfortunately, I still do not really understand. I have a step in the FV code where I update the state based on the fluxes. If you</div><div>violate CFD, this step is completely wrong.</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 dir="ltr"><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></blockquote><div><br></div><div>I can explain what I am doing, and maybe we can either change it to what you want, or allow you to use the pieces in your code to make things simpler.</div><div>Here is the comment from DMPlexComputeResidual_Internal():</div><div><br></div><div> /* FVM */<br> /* Get geometric data */<br> /* If using gradients */<br> /* Compute gradient data */<br> /* Loop over domain faces */<br> /* Count computational faces */<br> /* Reconstruct cell gradient */<br> /* Loop over domain cells */<br> /* Limit cell gradients */<br> /* Handle boundary values */<br> /* Loop over domain faces */<br> /* Read out field, centroid, normal, volume for each side of face */<br> /* Riemann solve over faces */<br> /* Loop over domain faces */<br> /* Accumulate fluxes to cells */<br></div></div></div></blockquote><div><br></div><div>Yes I saw the comments at the end of the method and they helped to understand.</div><div>Indeed at the end of this, you have the RHS of your system in the dU/dt = RHS sense. To be accurate it is the RHS at time t^n, like if it were a first order Euler approximation we would write U^(n+1) - U^n = dt * RHS^n.</div><div>Only with an implicit solver we want to write (U^(n+1) - U^n )/dt = RHS^(n+1) as you know. <br></div><div>Since we cannot really evaluate RHS^(n+1) we like to linearize it about RHS^n and we write RHS^(n+1) = RHS^n + [dRHS/dU]^n * (U^(n+1) - U^n)</div><div>All in all if we substitute we obtain something like [alpha * Identity - [dRHS/dU]^n](U^(n+1) - U^n) = RHS^n where alpha depends on the time step.</div><div>This left hand side matrix is what I would like to compute before advancing the system in time, whereas the RHS^n is what DMPlexComputeResidual_Internal gives.</div><div>Then there is the story of preconditioning this system to make the whole thing more stable but if I already could figure out how to generate that LHS it would be a great step in the right direction.</div><div><br></div>Thanks again a lot !</div><div class="gmail_quote"><br></div><div class="gmail_quote">Thibault<br></div><div class="gmail_quote"><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 dir="ltr"><div class="gmail_quote"><div></div><div>Obviously you might not need all the steps, and the last step is the one I cannot understand for your case.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</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>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><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>