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" <knep...@gmail.com>
À: "Ali ALI AHMAD" <ali.ali_ah...@utt.fr>
Cc: "petsc-maint" <petsc-ma...@mcs.anl.gov>, "petsc-users"
<petsc-users@mcs.anl.gov>, "Barry Smith" <bsm...@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_ah...@utt.fr |
ali.ali_ah...@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$
]
