<div dir="ltr"><div dir="ltr">On Thu, Dec 19, 2024 at 12:39 PM Aldo Bonfiglioli <<a href="mailto:aldo.bonfiglioli@unibas.it">aldo.bonfiglioli@unibas.it</a>> wrote:</div><div class="gmail_quote gmail_quote_container"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><u></u>
<div>
<p>Dear all,</p>
<p>I am using TS to solve the PDE </p>
<p>u_t + c * u_c = eps * u_xx in 1D i.e. smthg. similar to <a href="https://urldefense.us/v3/__https://petsc.org/release/src/ts/tutorials/ex4.c.html__;!!G_uCfscf7eWS!eeErx_FNg50N1eRlD9eyfoEp8I0llTQmVFpPWoz5_LQhL8K6uOId5zbgle0ghE71NJmon8T_vwwfN3LnEbf4cTUXtlkWJRflikU$" target="_blank">https://petsc.org/release/src/ts/tutorials/ex4.c.html</a></p>
<p>except that I am using DMPlex, instead of a structured grid.<br>
</p>
<p>Since the problem is linear, I tried smthg. like:</p>
<blockquote type="cite"><span style="font-family:monospace"><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
PetscCall(RHSMatrixHeat(ts,</span><span style="color:rgb(178,24,24);background-color:rgb(255,255,255)">0.0</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">,u,A,A,&appctx));
</span><br>
PetscCall(TSSetRHSFunction(ts,<span style="color:rgb(178,24,24);background-color:rgb(255,255,255)">NULL</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">,TSComputeRHSFunctionLinear,&appctx));
</span><br>
PetscCall(TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx));<br>
<br>
</span></blockquote>
<p>My problem is that the A matrix (as well as the global vector) is
missing the rows/cols corresponding to the endpoints where
Dirichlet boundary conditions are prescribed;</p>
<p> </p>
<blockquote type="cite">A global vector is missing both the shared
dofs which are not owned by this process, as well as <em>constrained</em>
dofs. These constraints represent essential (Dirichlet) boundary
conditions. They are dofs that have a given fixed value, so they
are present in local vectors for assembly purposes, but absent
from global vectors since they are never solved for during
algebraic solves.</blockquote>
<p>therefore, A*u = f(u) except in the two vertices adjacent to the
endpoints, i.e. first and last entries of the global vector.<br>
</p>
<p>The code is attached; when the flag -use_jacobian is set, f(u) is
computed as A*u, if it is not set, f(u) is computed in
FormRHSFunction.</p>
<p>The latter approach works, the former does not.<br>
</p>
<p></p></div></blockquote><div>I think I understand. When you declare Dirichlet conditions using Plex, it eliminates them from the solver. _However_, the</div><div>action must be computing using the local space, which contains the constraints. With FEM, this is natural. What you cannot</div><div>do is use an input vector in the global space, since it is missing the boundary values. Things you could possibly do:</div><div><br></div><div>1) Apply your operator in a matrix-free way, taking input in the full space, and giving output in the constrained space.</div><div> You can assemble the unconstrained operator if assembly is faster, and then map the output to the constrained space.</div><div><br></div><div>2) Ignore the constraint support in Plex and manage the constraints yourself using MatZeroRowsColumns() and constants</div><div> in the rhs for the boundary values. This is what DMDA tends to do.</div><div><br></div><div>3) Solve only for the update, since the boundary values do not change in the solve. This is what all the SNES/TS/TAO</div><div> support for linear solver using Plex does. We assemble using the local space (unconstrained), and output in the</div><div> global (constrained) space. This is why FormRHSFunction works.</div><div><br></div><div>Does that make sense?</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><p>Regards,</p>
<p>Aldo<br>
</p>
<pre cols="72">--
Dr. Aldo Bonfiglioli
Associate professor of Fluid Machines
Scuola di Ingegneria
Universita' della Basilicata
V.le dell'Ateneo lucano, 10 85100 Potenza ITALY
tel:+39.0971.205203 fax:+39.0971.205215
web: <a href="https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?m=002423__;!!G_uCfscf7eWS!eeErx_FNg50N1eRlD9eyfoEp8I0llTQmVFpPWoz5_LQhL8K6uOId5zbgle0ghE71NJmon8T_vwwfN3LnEbf4cTUXtlkWzhmq6Yg$" target="_blank">http://docenti.unibas.it/site/home/docente.html?m=002423</a></pre>
</div>
</blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><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="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bO7mL8fLGJ1HQeUpeYEkzSCNp7tF-J8HcWxslYKIr_sOi2FHpOJIG9c7fuoLoLEyWD7TSfN4tfhCDyC79CmM$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>