On Thu, Dec 19, 2024 at 12:39 PM Aldo Bonfiglioli < [email protected]> 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.... <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... <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_uCf... <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfs... >