[petsc-users] [petsc-maint] Question about NewtonTR / Inexact Newton with piecewise continuous functions
Barry Smith
bsmith at petsc.dev
Mon Aug 25 20:47:39 CDT 2025
You can put a debugger breakpoint in the SNESVI code that checks the bounds to see why it is not finding violations. This is done in SNESSolve_VINEWTONRSLS() which calls SNESVIGetActiveSetIS(). So check th variables in SNESVIGetActiveSetIS.
Note in even simple debuggers like gdb or lldb you can do
call VecView(x,0)
call VecView(xl,0)
etc
to examine the vectors quickly and verify if the values make sense.
Barry
> On Aug 25, 2025, at 6:02 PM, Ali ALI AHMAD <ali.ali_ahmad at utt.fr> wrote:
>
> Oh no, I was mistaken: I meant that the objective function is continuous, but its derivative is not. That is why I am using snesvinewtonssls.
>
> * However, the problem is that the nonlinear constraints I define in " ComputeBounds " are always inactive. As you can see in the output:
>
> 0 SNES Function norm 1.562978975444e+07
> 0 SNES VI Function norm 1.56298e+07 Active lower constraints 0/0 upper constraints 0/0 Percent of total 0. Percent of bounded 0.
> 0 KSP Residual norm 1.976350248286e+04
> 1 KSP Residual norm 1.483291315197e+04
> 2 KSP Residual norm 1.358944289069e+04
> 3 KSP Residual norm 1.276791956295e+04
> 4 KSP Residual norm 1.275078123378e+04
> 5 KSP Residual norm 1.253064296351e+04
> 6 KSP Residual norm 1.232954788760e+04
> 7 KSP Residual norm 1.232684700773e+04
> 8 KSP Residual norm 1.165265510468e+04
> 9 KSP Residual norm 1.159313964889e+04
> 10 KSP Residual norm 1.098241228864e+04
> 11 KSP Residual norm 1.063610434231e+04
> 12 KSP Residual norm 8.556981786904e+03
> 13 KSP Residual norm 7.792752698989e+03
> 14 KSP Residual norm 6.260039598596e+03
> 15 KSP Residual norm 3.025534878795e+03
> 16 KSP Residual norm 8.343990389000e+02
> 17 KSP Residual norm 2.591011944960e+02
> 18 KSP Residual norm 9.212497850148e+01
> 19 KSP Residual norm 4.137422005070e+01
> 20 KSP Residual norm 1.345189115925e+01
> 21 KSP Residual norm 2.795423327875e+00
> 22 KSP Residual norm 5.068536677416e-01
> 23 KSP Residual norm 1.167594367615e-01
> Linear solve converged due to CONVERGED_RTOL iterations 23
> [Non physique] pression négative ou NaN à la cellule 61 : p = -1.66778e+14
> Line search: objective function at lambdas = 1. is Inf or Nan, cutting lambda
> [Non physique] pression négative ou NaN à la cellule 61 : p = -4.06293e+13
> Line search: objective function at lambdas = 0.5 is Inf or Nan, cutting lambda
> [Non physique] pression négative ou NaN à la cellule 61 : p = -9.63498e+12
> Line search: objective function at lambdas = 0.25 is Inf or Nan, cutting lambda
> [Non physique] pression négative ou NaN à la cellule 61 : p = -2.15792e+12
> Line search: objective function at lambdas = 0.125 is Inf or Nan, cutting lambda
> [Non physique] pression négative ou NaN à la cellule 61 : p = -1568.11
> Line search: objective function at lambdas = 0.0625 is Inf or Nan, cutting lambda
> [Non physique] pression négative ou NaN à la cellule 788 : p = -11.2882
> Line search: objective function at lambdas = 0.03125 is Inf or Nan, cutting lambda
> Line search: Using full step: fnorm 1.562978975444e+07 gnorm 1.519904837982e+07
>
> So none of my lower or upper bounds are active during the iterations.
>
> * Here is my ComputeBounds function, where I try to enforce positivity of density and pressure:
>
> PetscErrorCode ComputeBounds(SNES snes, Vec xl, Vec xu) {
> PetscFunctionBeginUser;
>
> // 1) Get user context (must have been set via SNESSetApplicationContext)
> void *ctx = NULL;
> PetscCall(SNESGetApplicationContext(snes, &ctx));
> euler *e = (euler*)ctx;
>
> // Minimal positive thresholds
> const double eps_rho = 1e-8; // minimal density
> const double pmin = 1e-8; // minimal pressure
>
> // 2) Access current solution X (used here to build *dynamic* bounds)
> Vec X;
> const PetscScalar *x = NULL;
> PetscCall(SNESGetSolution(snes, &X));
> PetscCall(VecGetArrayRead(X, &x));
>
> // 3) Fill lower bounds 'xl'
> PetscScalar *l = NULL;
> PetscCall(VecGetArray(xl, &l));
>
> for (int K = 0; K < e->N_cells; ++K) {
> const int i_rho = Ncomp*K + 0; // density index
> const int i_m = Ncomp*K + 1; // momentum (1D). In 2D use mx,my and adjust kinetic term.
> const int i_rhoE = Ncomp*K + 2; // total energy index
>
> const double rho = (double)x[i_rho];
> const double m = (double)x[i_m];
>
> // Use a safe density to avoid division by ~0 in kinetic energy
> const double rho_safe = PetscMax(rho, eps_rho);
>
> // Kinetic energy per unit volume: |m|^2/(2*rho)
> const double kin_vol = 0.5 * (m*m) / rho_safe;
>
> // Enforce p >= pmin via total energy lower bound:
> // rhoE >= kinetic + p/(gamma-1)
> const double rhoEmin = kin_vol + pmin/(_gamma - 1.0);
>
> // Set lower bounds
> l[i_rho] = eps_rho; // rho >= eps_rho
> l[i_m] = PETSC_NINFINITY; // no lower bound on momentum
> l[i_rhoE] = rhoEmin; // rhoE >= rhoEmin (to keep p >= pmin)
> }
>
> // 4) Restore arrays
> PetscCall(VecRestoreArrayRead(X, &x));
> PetscCall(VecRestoreArray(xl, &l));
>
> // 5) Set upper bounds to +infinity (no upper constraints)
> PetscCall(VecSet(xu, PETSC_INFINITY));
>
> PetscFunctionReturn(PETSC_SUCCESS);
> }
>
> But when I run the solver, PETSc always reports Active lower constraints 0/0. Why are my constraints not activated?
>
> Best regards,
> ALI ALI AHMAD
>
> De: "Matthew Knepley" <knepley at gmail.com>
> À: "Ali ALI AHMAD" <ali.ali_ahmad at utt.fr>
> Cc: "petsc-maint" <petsc-maint at mcs.anl.gov>, "petsc-users" <petsc-users at mcs.anl.gov>, "Barry Smith" <bsmith at petsc.dev>
> Envoyé: Lundi 25 Août 2025 03:42:36
> Objet: Re: [petsc-maint] Question about NewtonTR / Inexact Newton with piecewise continuous functions
>
> On Sun, Aug 24, 2025 at 8:16 PM Ali ALI AHMAD <ali.ali_ahmad at utt.fr <mailto:ali.ali_ahmad at utt.fr>> wrote:
>> Hello,
>>
>> I am currently working with PETSc to solve a nonlinear system. My function is piecewise continuous (so it can be discontinuous), and its derivative (Jacobian) is also piecewise continuous.
>
> I am not aware of any convergence framework for piecewise continuous functions. They are not technically computable, meaning you could have to compute for a very very long time, and still not get an accurate output (since you have the jump). This also means that the derivative is not computable, and you can have arbitrarily large errors near the discontinuity. I don't see how you could prove convergence here, but maybe someone else knows something I don't.
>
> Thanks,
>
> Matt
>
>> I implemented the residual and also compute the Jacobian analytically.
>>
>> When the discontinuities are small, both Inexact Newton and NewtonTR (with scaled Newton direction) converge without problem.
>>
>> However, for test cases with larger discontinuities, sometimes the solver converges, but in other cases it fails to converge.
>>
>> I initially tried <vinewtonssls,vinewtonrsls> , but I cannot use them in my case, because my problem involves nonlinear constraints.
>>
>> So my question is: how does PETSc handle such situations internally (piecewise continuous objective/residual functions)?
>> And is there a recommended strategy within PETSc to deal with nonlinear solvers when and are discontinuous?
>>
>> I would like to continue working with PETSc and I am looking for a robust method to treat this type of problem.
>>
>> Thank you very much for your help and suggestions.
>>
>> Best regards,
>> ALI ALI AHMAD
>
>
>
> --
> 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!eVe1OReScpoGqh7jJoAWZ73himymGGPZO4ZjvamIdb0EmslI9gdViRqCNUS8J8Y1cimxSWK184MgrZJnsDic0eY$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!eVe1OReScpoGqh7jJoAWZ73himymGGPZO4ZjvamIdb0EmslI9gdViRqCNUS8J8Y1cimxSWK184MgrZJnYgB3K_I$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250825/213fa7a9/attachment-0001.html>
More information about the petsc-users
mailing list