<div dir="ltr"><div class="gmail_default" style="font-family:simsun,serif"><br></div><div class="gmail_default" style="font-family:simsun,serif">Thank you all! I think I did it. I used the ghost cells and the code behaves properly. Attached is my document with detailed derivation and the code.</div><div class="gmail_default" style="font-family:simsun,serif"><br></div><div class="gmail_default" style="font-family:simsun,serif">Best,</div><div class="gmail_default" style="font-family:simsun,serif">Qw</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Mark Adams <<a href="mailto:mfadams@lbl.gov">mfadams@lbl.gov</a>> 于2021年11月28日周日 下午10:17写道：<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">* I have used Matt's approach and it works but you do have to remember to update the ghost cells after you update your solution and before you use the operator again.<div>* I would use a convergence test to determine if the result is correct.</div><div>* Check that a constant function is a (trivial) solution to your discretization.</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sat, Nov 27, 2021 at 2:54 AM zhfreewill <<a href="mailto:zhfreewill@gmail.com" target="_blank">zhfreewill@gmail.com</a>> wrote:<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"><div class="gmail_default"><font face="tahoma, sans-serif">Hi,Barry, </font></div><div class="gmail_default"><font face="tahoma, sans-serif">Thanks for reply. </font><span style="font-family:tahoma,sans-serif">I meant I cannot get the right result as the FD version did, although the code did not report any error when running.</span></div><div class="gmail_default"><span style="font-family:tahoma,sans-serif"><br></span></div><div class="gmail_default"><font face="tahoma, sans-serif">Indeed, there were mistakes and I fixed it in the code (</font><font face="tahoma, sans-serif">I also attached additional derivation and code this time, please see the </font>attachment<span style="font-family:tahoma,sans-serif">). Still, I cannot get the right output so far. </span></div><div class="gmail_default"><br></div><div class="gmail_default"><font face="tahoma, sans-serif">Regards,</font></div><div class="gmail_default"><font face="tahoma, sans-serif">Qw</font></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> 于2021年11月26日周五 下午11:31写道：<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><div><br></div> Could you be more specific in saying how it fails (can you share the actual prototype code?) <div><br></div><div> Note here you seem to be missing a /dx term</div><div><br></div><div><blockquote style="font-family:Calibri,Helvetica,sans-serif;margin-top:0px;margin-bottom:0px"><div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = (u[i+1] - u[i]) - q_A; </b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153)"><b>//Neumann BC: du/dx(x=0.0) = q_A;</b></span></font></div><div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(153,153,153)"><b><br></b></span></font></div></blockquote><div><br></div><div> Barry</div><div><br><blockquote type="cite"><div>On Nov 26, 2021, at 4:02 AM, zhfreewill <<a href="mailto:zhfreewill@gmail.com" target="_blank">zhfreewill@gmail.com</a>> wrote:</div><br><div><div dir="ltr"><div style="font-family:simsun,serif"><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">Hi all,</span><br></div><div class="gmail_quote"><div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">I have a basic question on how to apply B.C. for a cell-centered structured FV discretization when using the SNES and Ts object.</span></div><div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt"><span class="gmail_default" style="font-family:simsun,serif"></span></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">For instance, for a 1-D transient diffusion equation</span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>u_t = u_xx, 0 < x < 1</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">with IC: </span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>u(x,t=0) = 0;</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">and Dirichlet BC:</span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>u(x=0) = u_A = 0.0;</b></span></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>u(x=1) = u_B = 1.0;</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">or Neumann BC:</span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>du/dx(x=0.0) = q_A = 0.0;</b></span></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>du/dx(x=1.0) = q_B = 1.0;</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">When using TS of PETSc, the equation can be organized in form of</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>U_t = G(t,u)</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">where
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>G</b></span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">
is the right<span class="gmail_default" style="font-family:simsun,serif"></span>-hand-side (RHS) function regardless of the specific method of discretization.</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Arial Black",Arial,sans-serif;font-size:10pt;color:rgb(222,106,25);background-color:rgba(0,0,0,0)"><b>(1) For FD discretization</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"><br>
</span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><b><span style="font-family:"Courier New",monospace;font-size:10pt"></span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)">A |---dx--| B</span></b></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>|-------|-------|-------|-------|-------| grid size = nx</b></span></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>0 i-1 i i+1 nx-1</b></span></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b><br>
</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">the vertex-centered FD discretization of the domain is</span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;font-size:10pt"><b></b></span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>x(i) = i*dx (i = 0,..., nx)), where dx = 1/nx</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">and the RHS function G of the equation is</span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;font-size:10pt"><b></b></span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>G[i] = (u[i+1]-2*u[i]+u[i-1])/dx^2</b></span></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b><br>
</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">There have been dozens of examples demonstrating way of BC setting with a finite difference medthod, e.g.,</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<a href="https://github.com/petsc/petsc/blob/386ebab930a846575609b49074a40c46e7f8ed75/src/ts/tutorials/ex26.c#L293-L343" id="gmail-m_-1618360744475031030gmail-m_-7161055579939979073gmail-m_1769095546239148101m_3788707207907009609LPlnk737449" target="_blank">https://github.com/petsc/petsc/blob/386ebab930a846575609b49074a40c46e7f8ed75/src/ts/tutorials/ex26.c#L293-L343</a><br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<div id="gmail-m_-1618360744475031030gmail-m_-7161055579939979073gmail-m_1769095546239148101m_3788707207907009609LPBorder_GTaHR0cHM6Ly9naXRodWIuY29tL3BldHNjL3BldHNjL2Jsb2IvMzg2ZWJhYjkzMGE4NDY1NzU2MDliNDkwNzRhNDBjNDZlN2Y4ZWQ3NS9zcmMvdHMvdHV0b3JpYWxzL2V4MjYuYyNMMjkzLUwzNDM." style="width:100%;margin-top:16px;margin-bottom:16px;max-width:800px;min-width:424px">
<div id="gmail-m_-1618360744475031030gmail-m_-7161055579939979073gmail-m_1769095546239148101m_3788707207907009609LPCloseButtonContainer950726" title="删除链接预览" role="button">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">which I think I hav</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">e underst</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">ood
</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">how it works and the following lines of code evaluate the RHS function
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(200,38,19);background-color:rgba(0,0,0,0)"><b>G(t,u)</b></span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">
in PETSc's </span><b><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(200,38,19);background-color:rgba(0,0,0,0)">TSSetRHSFun</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(200,38,19);background-color:rgba(0,0,0,0)">ction()</span></b><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">
routine</span></div>
</div>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif">
<font size="1"><span style="font-family:"Courier New",monospace"></span><span style="font-family:"Courier New",monospace;color:rgb(237,92,87);background-color:rgba(0,0,0,0)">/*------------------FD discretization-------------------------*/</span></font></div>
<div style="font-family:Calibri,Helvetica,sans-serif">
<span style="font-family:"Courier New",monospace;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b><font size="1">for (i = xs; i < xs + xm; i++) {</font></b></span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;color:rgb(153,153,153)"><b><font size="1">// left boundary</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">if (i == 0) {</font></b></span></div>
<blockquote style="margin-top:0px;margin-bottom:0px">
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = u[i] - u_A; </b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//
Dirichlet BC: u(x=0) = u_A;</b></span></font></div>
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = (u[i] - u[i+1])/dx - q_A; </b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//
Neumann BC: du/dx(x=0.0) = q_A;</b></span></font></div>
</blockquote>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">}
</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b><font size="1">// right boundary</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">else if (i == mx - 1) {</font></b></span></div>
<blockquote style="margin-top:0px;margin-bottom:0px">
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = u[i] - u_B; </b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//
Dirichlet BC: u(x=1.0) = u_B;</b></span></font></div>
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = (u[i]-u[i-1])/dx - q_B; </b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//
Neumann BC: du/dx(x=0.0) = q_A;</b></span></font></div>
</blockquote>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">}</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b><font size="1">// internal</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">else {</font></b></span></div>
<blockquote style="margin-top:0px;margin-bottom:0px">
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = (u[i-1] - 2 * u[i] + u[i + 1])/dx^2;
</b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//u_xx</b></span></font></div>
</blockquote>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">}</font></b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif">
<font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>}</b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><span style="color:rgb(237,92,87);font-variant-ligatures:inherit;font-variant-caps:inherit"><b>/*------------------------------------------------------------*/</b></span></span></font></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;color:rgb(237,92,87);font-size:8pt"><b><br>
</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">The BCs, e.g., at the left boudnary, can be straightforwardly set as follows:</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">Dirichlet BC:</span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>G[0] = u[0] - u_A</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">which indicates</span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>G[0] = u[0] - u_A = 0</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">and Dirichlet BC</span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>u[0] = u_A</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">is thus applied</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">Neumann BC:</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>G[0] = (u[0] - u[1])/dx - q_A</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">which indicates</span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>G[0] = (u[0] - u[1])/dx - q_A = 0</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">and Neumann BC</span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>du/dx = (u[0] - u[1])/dx = q_A</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">is applied</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">This is the way the above examples apply BCs.</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">For a structured, cell-centered FVM discretization, however, I tried and found the BCs can not be applied in a similar way. I will explain how I tried this as follows.</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Arial Black",Arial,sans-serif;font-size:10pt;color:rgb(222,106,25)"><b>(2) FV discretization</b></span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><b><span style="font-family:"Courier New",monospace;font-size:10pt"></span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)">A |---dx--| B</span></b></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>|---o---|---o---|---o---|---o---|---o---| cell size = nx</b></span></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>W w P e E
</b></span></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>0 i-1 i i+1 nx-1</b></span></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b><br>
</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">(NOTE: the center cell is denoted as P and its neighbouring cells as W, E; </span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">cell faces are denoted as w, e; the left
and right boundary faces are A and B, respectively.)</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">The cell-centered FV discretization of the domain is</span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;font-size:10pt"><b></b></span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>x(i) = dx/2 + i*dx (i = 0,..., nx-1)), where dx = 1/nx</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">the diffusion term u_xx can be reorganized by using divergence theorem</span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><b><span style="font-family:"Courier New",monospace;font-size:10pt"></span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)">(u_E-u_P) (u_P-u_W)</span></b></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>--------- - ---------</b></span></div>
<div><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b> dx dx</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">and the RHS function G of the equation is</span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;font-size:10pt"><b></b></span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>G[i] = (u_E-u_P)/dx - (u_P-u_W)/dx = (u[i+1]-2*u[i]+u[i-1])/dx</b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"><span style="margin:0px;font-size:10pt;background-color:rgb(255,255,255);font-family:Arial,Helvetica,sans-serif">the key lines of code in RHS function G(t,u) in PETSc's</span><span style="margin:0px;font-size:10pt;background-color:rgb(255,255,255)"><span style="font-family:Arial,Helvetica,sans-serif"> </span></span><span style="margin:0px;font-size:10pt;color:rgb(237,92,87);background-color:rgb(255,255,255);font-family:Arial,Helvetica,sans-serif"><b>TSSetRHSFunction()</b></span><span style="margin:0px;font-size:10pt;background-color:rgb(255,255,255)"><span style="font-family:Arial,Helvetica,sans-serif"> </span></span><span style="margin:0px;font-size:10pt;background-color:rgb(255,255,255);font-family:Arial,Helvetica,sans-serif">routine
can be</span></span></div><div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"><span style="font-family:"Courier New",monospace;font-size:10pt"><span style="margin:0px;font-size:10pt;background-color:rgb(255,255,255);font-family:Arial,Helvetica,sans-serif"><br></span></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif">
<b><font size="1"><span style="font-family:"Courier New",monospace"></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)">/*------------------FV discretization-------------------------*/</span></font></b></div>
<div style="font-family:Calibri,Helvetica,sans-serif">
<span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">for (i = xs; i < xs + xm; i++) {</font></b></span></div>
<blockquote style="font-family:Calibri,Helvetica,sans-serif;margin-top:0px;margin-bottom:0px">
<div><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b><font size="1">// left boundary</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">if (i == 0) {</font></b></span></div>
<blockquote style="margin-top:0px;margin-bottom:0px">
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = 3.0 * u[i] - u[i+1] - 2.0 * u_A;
</b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//Dirichlet BC: u(x=0.0) = u_A, does not work</b></span></font></div>
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = (u[i+1] - u[i]) - q_A;
</b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//Neumann BC: du/dx(x=0.0) = q_A;</b></span></font></div>
</blockquote>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">}</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b><font size="1">// right boundary</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">else if (i == mx - 1) {</font></b></span></div>
<blockquote style="margin-top:0px;margin-bottom:0px">
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = 3.0 * u[i] - u[i-1] - 2.0 * u_B;
</b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//Dirichlet BC: u(x=1.0) = u_B, does not work</b></span></font></div>
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = (u[i]-u[i-1])/dx - q_B; </b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//Neumann
BC: du/dx(x=1.0) = q_B, does not work</b></span></font></div>
</blockquote>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">}</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b><font size="1">// internal</font></b></span></div>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">else {</font></b></span></div>
<blockquote style="margin-top:0px;margin-bottom:0px">
<div><font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>G[i] = (u[i-1] - 2*u[i] + u[i+1])/dx; </b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><b>//u_x</b></span></font></div>
</blockquote>
<div><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b><font size="1">}</font></b></span></div>
</blockquote>
<div style="font-family:Calibri,Helvetica,sans-serif">
<font size="1"><span style="font-family:"Courier New",monospace;color:rgb(237,92,87)"><b>}</b></span><span style="font-family:"Courier New",monospace;color:rgb(153,153,153);background-color:rgba(0,0,0,0)"><span style="color:rgb(237,92,87);font-variant-ligatures:inherit;font-variant-caps:inherit"><b>/*------------------------------------------------------------*/</b></span></span></font></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">I tried to apply BCs in a way similar to that of the FD discretization.</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">For the Dirichlet BC on the left boundary, because of the cell-centered arrangement of u[i], the left BC is given as
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>u_A</b></span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt"> at cell face A instead of cell center.</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">So, I first applied a linear interpolation to get the cell-centered value of
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>u[0]</b></span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt"> from the left boundary value
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>u_A</b></span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt"> and neighbouring
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>u[1]</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>u[0] = (u[1]+2*u_A)/3.0 (equivalent to 3*u[0] - u[1] - 2*u_A = 0.0)</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b> |-dx/2-|------dx-----|</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b> |------o------|------o------| linear interpolation: u[0] = (u[1]+2*u_A)/3.0</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b> A P e E</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b> u_A u[0] u[1]</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">Dirichlet BC:</span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>G[0] = 3.0*u[0] - u[1] - 2.0*u_A</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> in which </span>
<span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>u[0] = (u[1]+2*u_A)/3.0</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">is implicitly calculated and I expect the Dirichlet BC</span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>u(x=0) = u_A</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">is applied</span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">For the Neumann BC, left BC is given as a gradient of
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>u</b></span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">,
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87);background-color:rgba(0,0,0,0)"><b>q_A</b></span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">, at cell face A</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">Neumann BC:</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>G[i] = (u[1] - u[0])/dx - q_A</b></span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:"Courier New",monospace;font-size:10pt"> </span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">in which Neumann BC</span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:"Courier New",monospace;font-size:10pt;color:rgb(237,92,87)"><b>du/dx = (u[1] - u[0])/dx = q_A</b></span><span style="font-family:"Courier New",monospace;font-size:10pt">
</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt;background-color:rgba(0,0,0,0)">is expected to be applied.</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">However, this FV discretization failed to work. Is there anything wrong or c</span><span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">ould you give me a hint?</span></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">Best,</span></div>
<span style="font-family:Arial,Helvetica,sans-serif;font-size:10pt">Qw</span><br>
</div>
</div>
</div>
</div>
</div>
</div>
</div></div>
</div></blockquote></div><br></div></div></blockquote></div></div>
</blockquote></div>
</blockquote></div>