[petsc-users] [petsc-maint] Question about NewtonTR / Inexact Newton with piecewise continuous functions

Mon Aug 25 17:02:12 CDT 2025

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 < [ mailto:ali.ali_ahmad at utt.fr | 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 

BQ_BEGIN

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 

BQ_END

-- 
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/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!c59r-_cWhRyXgrv6zPsyXOKKCkQh6vxCZENKwjHfU9j8W-SbCNE9hvNJQKo-F6mc6MAmkFUmr5Ysk_S5VvLlpZfDMhUtWg$  | https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!c59r-_cWhRyXgrv6zPsyXOKKCkQh6vxCZENKwjHfU9j8W-SbCNE9hvNJQKo-F6mc6MAmkFUmr5Ysk_S5VvLlpZeH41xL_Q$  ] 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250826/c4f79c51/attachment.html>