[petsc-users] Advice on setting Dirichlet bcs in DMPlex

[Matthew Knepley]

On Thu, Dec 19, 2024 at 12:39 PM Aldo Bonfiglioli <
aldo.bonfiglioli at unibas.it> wrote:

> Dear all,
>
> I am using TS to solve the PDE
>
> u_t + c * u_c = eps * u_xx in 1D i.e. smthg. similar to
> https://urldefense.us/v3/__https://petsc.org/release/src/ts/tutorials/ex4.c.html__;!!G_uCfscf7eWS!bO7mL8fLGJ1HQeUpeYEkzSCNp7tF-J8HcWxslYKIr_sOi2FHpOJIG9c7fuoLoLEyWD7TSfN4tfhCD1mcIe-Z$ 
> <https://urldefense.us/v3/__https://petsc.org/release/src/ts/tutorials/ex4.c.html__;!!G_uCfscf7eWS!eeErx_FNg50N1eRlD9eyfoEp8I0llTQmVFpPWoz5_LQhL8K6uOId5zbgle0ghE71NJmon8T_vwwfN3LnEbf4cTUXtlkWJRflikU$>
>
> except that I am using DMPlex, instead of a structured grid.
>
> Since the problem is linear, I tried smthg. like:
>
>    PetscCall(RHSMatrixHeat(ts,0.0,u,A,A,&appctx));
>    PetscCall(TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx));
>
>
>    PetscCall(TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx));
>
> 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;
>
> A global vector is missing both the shared dofs which are not owned by
> this process, as well as *constrained* 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.
>
> therefore, A*u = f(u) except in the two vertices adjacent to the
> endpoints, i.e. first and last entries of the global vector.
>
> 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.
>
> The latter approach works, the former does not.
>
> I think I understand. When you declare Dirichlet conditions using Plex, it
eliminates them from the solver. _However_, the
action must be computing using the local space, which contains the
constraints. With FEM, this is natural. What you cannot
do is use an input vector in the global space, since it is missing the
boundary values. Things you could possibly do:

1) Apply your operator in a matrix-free way, taking input in the full
space, and giving output in the constrained space.
    You can assemble the unconstrained operator if assembly is faster, and
then map the output to the constrained space.

2) Ignore the constraint support in Plex and manage the constraints
yourself using MatZeroRowsColumns() and constants
     in the rhs for the boundary values. This is what DMDA tends to do.

3) Solve only for the update, since the boundary values do not change in
the solve. This is what all the SNES/TS/TAO
    support for linear solver using Plex does. We assemble using the local
space (unconstrained), and output in the
    global (constrained) space. This is why FormRHSFunction works.

Does that make sense?

  Thanks,

     Matt

> Regards,
>
> Aldo
>
> --
> 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: https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?m=002423__;!!G_uCfscf7eWS!bO7mL8fLGJ1HQeUpeYEkzSCNp7tF-J8HcWxslYKIr_sOi2FHpOJIG9c7fuoLoLEyWD7TSfN4tfhCD1nxZp-t$  <https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?m=002423__;!!G_uCfscf7eWS!eeErx_FNg50N1eRlD9eyfoEp8I0llTQmVFpPWoz5_LQhL8K6uOId5zbgle0ghE71NJmon8T_vwwfN3LnEbf4cTUXtlkWzhmq6Yg$>
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bO7mL8fLGJ1HQeUpeYEkzSCNp7tF-J8HcWxslYKIr_sOi2FHpOJIG9c7fuoLoLEyWD7TSfN4tfhCD4BXqfL8$  <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bO7mL8fLGJ1HQeUpeYEkzSCNp7tF-J8HcWxslYKIr_sOi2FHpOJIG9c7fuoLoLEyWD7TSfN4tfhCDyC79CmM$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20241219/9a2c3017/attachment.html>