<div dir="ltr"><div class="gmail_default" style="font-family:simsun,serif">Thank you for your reply. This is a good FD example for me and I think I understood it at least for FD discretization.<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">Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> 于2021年11月16日周二 下午7:11写道：<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 Mon, Nov 15, 2021 at 10:24 PM zhfreewill <<a href="mailto:zhfreewill@gmail.com" target="_blank">zhfreewill@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 style="font-family:simsun,serif">Hi, <br>I'm new to PETSc. I have been confused about how to apply boundary conditions (BCs) (e.g., the Dirichlet and the Neumann type) when using the time-stepping (TS) object of PETSc.</div><div style="font-family:simsun,serif"><br>I have read the example codes (<a href="https://petsc.org/release/src/ts/tutorials/index.html" target="_blank">https://petsc.org/release/src/ts/tutorials/index.html</a>) and vaguely found that BCs seems to be set as follows:<br><br>1. when evaluating a matrix of the right hand side function g(x) (e.g., using a function like RHSMatrixHeat)<br>2. when evaluating a matrix of the Jacobian matrix of g'(x) (e.g., using a function like FormJacobian)<br></div><div style="font-family:simsun,serif"><br>For instance, in the ex1.4 (<a href="https://petsc.org/release/src/ts/tutorials/ex4.c.html" target="_blank">https://petsc.org/release/src/ts/tutorials/ex4.c.html</a>), with comments at the beginning states<br><br> 17: /* ------------------------------------------------------------------------<br> 19: This program solves the one-dimensional heat equation (also called the<br> 20: diffusion equation),<br> 21: u_t = u_xx,<br> 22: on the domain 0 <= x <= 1, with the boundary conditions<br> 23: u(t,0) = 0, u(t,1) = 0,<br> 24: and the initial condition<br> 25: u(0,x) = sin(6*pi*x) + 3*sin(2*pi*x).<br> 26: This is a linear, second-order, parabolic equation.<br><br> 28: We discretize the right-hand side using finite differences with<br> 29: uniform grid spacing h:<br> 30: u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2)<br>...<br> 47: ------------------------------------------------------------------------- */<br><br> and Dirichlet BCs, u(t,0) = 0 and u(t,1) = 0, are to set on the both side of the discretized domain, the corresponding codes can be found on the following lines in a function RHSMatrixHeat:<br> <br>/* ------------------------------------------------------------------------- */</div><div style="font-family:simsun,serif">463: RHSMatrixHeat - User-provided routine to compute the right-hand-side matrix for the heat equation. */<br>491: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,void *ctx)<br>492: {<br>...<br>505: /* Set matrix rows corresponding to boundary data */<br>509: if (mstart == 0) { /* first processor only */<br>510: v[0] = 1.0;<br>511: MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);<br>512: mstart++;<br>513: }<br>515: if (mend == appctx->m) { /* last processor only */<br>516: mend--;<br>517: v[0] = 1.0;<br>518: MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);<br>519: }<br><br>521: /* Set matrix rows corresponding to interior data. We construct the matrix one row at a time. */<br>525: v[0] = sone; v[1] = stwo; v[2] = sone;<br>526: for (i=mstart; i<mend; i++) {<br>527: idx[0] = i-1; idx[1] = i; idx[2] = i+1;<br>528: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);<br>529: }<br>...<br>550: }</div><div style="font-family:simsun,serif">/* ------------------------------------------------------------------------- */</div><div style="font-family:simsun,serif"><br><br>My questions are:<br><br>1. How do these lines of code embody the Dirichlet BCs u(t,0) = 0 and u(t,1) = 0.<br></div></div></blockquote><div><br></div><div>This is the finite difference version of BC. It is just moving the term for a known value to the forcing. However, since</div><div>the BC is 0, then we just do not compute the term. If the BC were nonzero, you would also see the term in the</div><div>residual. For example,</div><div><br></div><div> <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex19.c#L271" target="_blank">https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex19.c#L271</a></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 style="font-family:simsun,serif">2. Similarly, if I want to apply a Neumann type BCs in these lines, How can I do?<br></div></div></blockquote><div><br></div><div>For finite differences, you usually just approximate the derivative as the one-sided difference of the last two values,</div><div>and set this equal to the known value. For finite elements, Neumann conditions are just an addition weak form</div><div>integrated over the boundary.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div style="font-family:simsun,serif">Any suggestions or documents?<br><br>Best,<br>Qw<br></div></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>