Repository navigation
Conversation
Largely had Claude drive this port as it the original SNLS trust region dogleg solver was relatively straightforward. Most of it would map over to MFEM's framework and wasn't overly complicated. The big boon though from having Claude do this is I also had it add the necessary mult transpose operators for the nonlinearform integrators and in particular this was done for the PA forms as the EA was trivial and the full matrix version already implemented it. I was also surprised that it was able to do the BBar PA Grad implementation of everything as well. The math behind it was actually more straightforward than I was expecting, and it was cool to see how it could derive all of it using different relationship it had discovered from looking at the full integration form.
… mfem::Array<bool> Was getting some odd GPU failures here at one point and just moved over to std::array<bool, 3> to get rid of them as I didn't want to deal with the odd memory issues I was hitting. The MFEM stuff was still fine. I just think I was hitting an odd bug in some GPU kernel, but moving to the std::array is ultimately the better approach as it honestly makes more sense for this stuff...
Had codex help work out what changes were needed to get ExaConstit updated to work with MFEM v4.9+ as I wanted to make sure we could use newer versions of Hypre as apparently that works better with ROCm v7+
Found a couple bugs in the post-processing: First one is that volume average values were getting outputted every time step designated for the viz files and the viz files were getting outputted every time step. So caught and fixed that issue. Next found in certain cases there could be segfaults due to some models not having a variable defined and thus a vdim equal to 0 causing segfaults which was a fun bug to run down... Finally this was one I had codex help dissect and chase down, we had some fun MPI stalls occurring and it was due to how data was being solved requiring the global communicator rather than a region defined communicator... Didn't notice this in earlier testing as our regions were usually on every rank... I'd still blame this one though on how MFEM is handling the communicators for meshes / par finite element spaces and in turn how this relates to the data collections as well...
Was running into a nasty GPU-bug where on newer versions of MFEM if we used > 1 GPUs we got different answers than our CPU runs as certain terms were just 0.0. I could not for the life of me figure it out other than it was likely due to some thing in the velocity field near the time the boundary conditions were being applied was not getting set. I threw codex at the problem and it was able to over a couple iterations of debugging work out the error and find a suitable new MFEM API that we could use that fixed the GPU error and kept our answers on the CPU the same.
There were a few regressions in handling orientation and grain files in the options as well as where to put the orientations with respect to ECMech models. Had codex help fix this issue as I've been busy trying to get other features worked out that are taking a bit higher priority. Although, the solution proposed was simple enough and not terribly far off from what I might have done.
…t following what it should have
…andled in more recent version of clang-format
…, plus TRDOG rejection fixes
The trust-region dogleg (TRDOG) solver printed its progress in its own style
("TR dogleg: iter N, ||r|| = ..., delta = ... [DL]") along with a scattering of
"TR: rho = ...", "TR: increased delta to ...", and "TR dogleg: stepping along first
leg ..." lines from inside TrDeltaControl and Dogleg(), none of which lined up with what
the NR and NRLS solvers print. Separately, a user asked for a way to quiet the
per-iteration nonlinear solver output. Both are addressed here by giving all three
solvers one shared output path and a single verbosity option, so a single parser can
follow how any of our nonlinear solvers are doing over a simulation.
New `print_info` integer in [Solvers.NR] (NonlinearSolverOptions::print_info, parsed in
option_solvers.cpp, validated to be >= 0, shown in the options printout, documented in
options.toml). It is passed straight to SetPrintLevel(), replacing the hard-coded
SetPrintLevel(1) in SystemDriver. That keeps the existing rank-0-only behavior since
MFEM's IterativeSolver::SetPrintLevel sets print_level to -1 on every other rank.
- 0: only a one line summary once the solver converges or fails, giving the number of
iterations, the final residual, and the final relative residual. Some basic info is
always printed rather than none at all.
- 1: residual and step scale every iteration (default; NR/NRLS output is the same
information as before, just in the unified format below).
- 2+: also prints a detailed per-trial line for the trust-region solver (same as 1 for
NR and NRLS).
Three protected helpers on ExaNewtonSolver are used by all three solvers:
PrintIterationInfo, PrintStepScaleInfo, and PrintSummaryInfo. Every per-iteration line
now has the form "Newton iteration N : key = value, ...":
Newton iteration 0 : ||r|| = 3.62889e-05
Newton iteration 0 : step scale = 0.0272538
Newton iteration 1 : ||r|| = 3.0976e-05, ||r||/||r_0|| = 0.853593
and the level 0 summary is:
Newton solver converged in 14 iterations : ||r|| = 2.04178e-07, ||r||/||r_0|| = 5.35934e-06
The step scale on iteration N is the fraction of the full Newton step taken from that
iteration, and each solver supplies its own version of it: the 1.0/0.5 relaxation
factor for NR (now printed every iteration rather than only when it was reduced), the
line search factor for NRLS (replacing "The relaxation factor for this iteration is"),
and for TRDOG the fraction of the Newton step the trust region allows (1 for a full
Newton step, delta / ||s_N|| otherwise, since a non-Newton dogleg step lies on the
trust-region boundary). A step scale that gets small and stays small is an early sign
of trouble no matter which solver is running.
- TrDeltaControl (DecrDelta / IncrDelta / UpdateDelta) and Dogleg() no longer print
anything; the print_level parameters on the TrDeltaControl methods are removed.
Dogleg() now reports which leg it took through a new DoglegStep enum (NEWTON,
STEEPEST_DESCENT, DOGLEG) along with beta, replacing the use_nr out parameter.
- Only accepted steps count as iterations now. A rejected trial is retried from the same
state with a smaller trust region, reusing the Jacobian and Newton step, so it costs
one residual evaluation and no Krylov solve. That is the same kind of work as NRLS's
extra line-search evaluations, which have never counted as iterations. This keeps
max_iter, final_iter (and therefore the auto time step's dt control), and the
"Newton iteration N" lines meaning the same thing for all three solvers. This is a
deliberate difference from SNLS, where rejected steps count toward _nIters; the
trust-region path itself (trial steps, delta updates, accept/reject decisions) is
unchanged.
- At print_info = 1, rejected trials print nothing, the same way NRLS's line-search
trials print nothing; only the accepted step's scale and resulting residual are
shown.
- At print_info >= 2, each trial (accepted or rejected) prints a single indented line
under the iteration it was taken from:
TR step: dogleg (beta = 0.172956), predicted ||r|| = 0.0039825, trial ||r|| = 0.0035787, rho = 1.24219, delta = 0.00086306 -> 0.00129459, accepted
These lines do not start with "Newton iteration", so a parser of the unified format
can simply skip them.
- The level 0 summary gains ", rejected steps = K" for TRDOG so the extra work done is
still visible when the per-iteration output is off.
- A hidden internal cap of 50 consecutive rejections stops the solve as a safeguard
against an endless loop. It is not a user option; in practice deltaMin is reached
well before it.
- Fixed in passing: when a rejected trial also hit a delta-control failure, the loop
stopped before restoring x, so the solver reported the rejected trial's larger
residual and left x at that point. The previous state is now restored first.
After a rejected full Newton step, DecrDelta's took_full branch sets
delta = xiDecDelta * sqrt(delta * ||s_N||). The full Newton step was taken because
||s_N|| <= delta, and the new delta is still >= ||s_N|| whenever
delta >= ||s_N|| / xiDecDelta^2 (16 * ||s_N|| with the default xiDecDelta = 0.25). In that
case the retry, made from the same restored x with the same Newton step, is the exact
same step: same residual, same rho, and the same rejection again, wasting a full
residual evaluation (a constitutive update at every quadrature point) each time.
Visible directly in the new level 2 output:
TR step: Newton, ..., trial ||r|| = 0.0118788, rho = -0.63268, delta = 1 -> 0.022842, rejected
TR step: Newton, ..., trial ||r|| = 0.0118788, rho = -0.63268, delta = 0.022842 -> 0.00345224, rejected
UpdateDelta now keeps calling DecrDelta after a rejected full Newton step until
delta < ||s_N|| (or deltaMin is reached). The skipped retries were deterministic
repeats, so this lands on exactly the same delta the wasted evaluations would have
reached, and the solve follows the same path with fewer residual evaluations. The loop
always terminates: DecrDelta strictly decreases delta while delta >= ||s_N||, and
repeating it converges to xiDecDelta^2 * ||s_N||. SNLS's own TrDeltaControl::decrDelta
has the same issue (both updateDelta variants, plus the rjSuccess == false paths in
SNLS_kernels.h and SNLS_kernels_batch.h) and will be fixed there separately.
Build is clean, and the full src tree was run through .clang-format per
developers_guide.md. All runs below are on 2 MPI ranks, confirming only rank 0 prints.
Ran voce_full.toml-based cases (scratch copies, with the result folders removed
afterwards) at every print_info level for NR, NRLS, and TRDOG. The normal problem only
ever takes full Newton steps, so a tiny delta_init was forced to exercise the steepest
descent and dogleg legs, and a 1-iteration limit was used to check the failure summary.
Rejections were exercised with fixed dt = 5 and dt = 20 single-step runs. For dt = 5 the
accepted-iteration output (14 iterations, every residual and step scale) is
line-for-line identical before and after the DecrDelta fix, while rejected steps drop
from 4 to 3 (6 to 5 for dt = 20) -- the duplicate trial is gone and nothing else
changed. ctest: the same pre-existing failures as before these changes and nothing
else -- the multi_material_test / voce_full_cyclic* near-zero noise differences.
…a symlink On HPC systems the results directory is commonly a symlink to another file system (`ln -s /remote/dir/location results`), and runs there occasionally ended up with an odd directory layout, e.g. a real `results/<basename>` alongside a `results/test_case -> ~/test_case` link. The output path was being resolved over and over: every rank, every PostProcessingFileManager instance, the logger, and each file write called weakly_canonical() on it independently. On top of that, the code paired checks with later actions (exists() then create_directories(), exists() to decide whether to write a header, etc.). With a symlink in the path, and especially with several ExaConstit runs on the same node creating or changing that symlink, different ranks and different points in a run could see different states of the path between a check and the action that relied on it (a TOCTOU race). libstdc++/libc++ do not make these filesystem calls atomic with respect to each other, so this has to be handled on our side. ExaOptions::parse_options now ends with a one-time collective step for post_processing.volume_averages.output_directory: - Rank 0 makes the path absolute, creates it with create_directories() (which follows a symlink already in place, or otherwise creates a real directory), then resolves it with canonical() to the physical directory that now exists. - The resolved path is broadcast to all ranks and replaces the option value, so every later consumer (post-processing, visualization data collections, logger) works from the same absolute physical path. - Creating before resolving means nothing after this point goes back through a symlink in the user-provided path. Another run creating, retargeting or removing the `results` link later on cannot redirect where this run writes. If the path changes in the small window between the two calls, canonical() fails rather than silently resolving somewhere unexpected. - All earlier error paths in parse_options abort on every rank, so every rank always reaches the broadcast. Following from that: - PostProcessingFileManager's constructor, EnsureDirectoryExists, CreateOutputFile, and UnifiedLogger::initialize / end_capture use the path as given; every weakly_canonical() call is gone. - EnsureDirectoryExists no longer re-resolves and re-broadcasts the path. Rank 0 calls create_directories() directly (it is a no-op for an existing directory), checks is_directory(), does the existing write test, and broadcasts success. Its message still says "Creating output directory" or "Using existing directory" depending on what create_directories() reported. - The logger's log directory and captured-output directories are likewise created without an exists() check first. - WriteVolumeAverage decides whether to write the header by checking whether the opened file is empty (seekp to end, tellp() == 0) instead of calling exists() before opening it. A results link created before its target directory exists (`ln -s /scratch/new_dir results`) is a realistic setup, but create_directories() does not create a dangling link's target -- it just fails because the link itself exists. Previously that meant the run went through the whole simulation, exited with 0, and saved nothing, with only a series of warnings to show for it. A new create_dangling_symlink_targets() helper walks the output path and creates the target of every dangling link along it (relative targets, deeper missing paths, and links that point through other links, with a depth cap guarding against cycles). If the output directory still can't be set up after that (for example `results` is a regular file, or a link to one), every rank now aborts right away with an error naming the directory, instead of running the whole simulation with nowhere to write. CreateOutputFile called MPI_Barrier(comm) whenever the target directory looked missing. It is only called by the writing rank (WriteVolumeAverage on rank 0 of the region communicator), so the other ranks never reached that barrier. It has been removed; the directory creation does not need any synchronization there. Build is clean, and the full src tree was run through .clang-format per developers_guide.md. Reproduced the user setup (`results -> ~/test_case`) with 8-rank voce_ea runs. Output lands in the link target and no extra links or nested directories appear. A run whose link was retargeted mid-run kept writing all of its output to the original target while a second run started afterwards wrote to the new one. Stress tested with a script that runs 2-6 concurrent simulations (12 ranks total, 6/4/3/2 ranks each) from the same directory, all writing to a relative `results`, starting from each of: missing, real directory, symlink, dangling symlink, and missing with a link or real directory created at a random moment after the runs start. In about half the trials a background thread also kept changing `results` every 2-80 ms while the runs were going (atomic retargets, including to targets that don't exist yet, unlinks, new links, new real directories). Each run was checked for: - a normal exit within the time limit; - all output in exactly one physical location; - every average file having 1 header and the full set of rows; - every VisIt cycle being present; - stresses matching a reference run with the same rank count to 1e-12. Each trial was also checked for stray links or nesting. Two rounds of 150 trials (1208 runs, 107 of them through dangling links) all passed with no hangs. Dangling-link cases (absolute, deep, relative, chained, and concurrent) all create their targets and run normally. `results` being a regular file, or a link to one, aborts cleanly on all ranks in under a second. ctest: no new failures from this change, and the full suite passes in the CI environment (Ubuntu 24.04, GCC 13, ExaCMech v0.4.3) on both x86_64 and aarch64.
The two Rust/PyO3 Python extensions in the repo (xtal_light_up and voxel_coarsen) were
both on pyo3 0.21, which has two RustSec advisories that Dependabot flags:
- RUSTSEC-2025-0020: risk of buffer overflow in PyString::from_object (fixed in 0.24.1)
- RUSTSEC-2026-0177: missing Sync bound on PyCFunction::new_closure closures (fixed in
0.29.0)
Both crates now use pyo3 0.29 and the matching numpy 0.29, plus anyhow 1.0.104 (1.0.86
had an unsoundness warning in Error::downcast_mut, RUSTSEC-2026-0190).
- pyo3 / numpy 0.21 -> 0.29, anyhow -> 1.0.104, Cargo.lock refreshed. cargo audit is
now clean.
- #[pyfn] is deprecated in 0.29 and is going away in a future PyO3 release, so the seven
functions are now top-level #[pyfunction]s registered in the #[pymodule] with
wrap_pyfunction!. Python names and signatures are unchanged. Most of the diff in
mod.rs is the 4-space dedent from moving the functions out of the module function;
their bodies are untouched.
- pyo3 / numpy 0.21 -> 0.29, anyhow -> 1.0.104.
- abi3-py37 -> abi3-py39 since PyO3 0.29 no longer supports Python 3.7 (matches
xtal_light_up).
- The optional polars dependency goes 0.30 -> 0.55. 0.30 pulled in arrow2, fast-float
and lexical-core advisories with no fixed version on that line, and it no longer
compiled with current Rust anyway.
Fixed in passing, as both optional features were broken:
- python feature: did not compile due to a `Bound<\`py` typo in pycoarsen/mod.rs. The
binding also moves from #[pyfn] to #[pyfunction] + wrap_pyfunction!, and
from_vec_bound -> from_vec for the numpy 0.29 API.
- polars feature: read_data moves to the current CsvReadOptions API, and its delimiter
and comment characters were swapped ('%' as the delimiter, ',' as the comment). It
now uses ',' and '%' like the default data_reader path.
- README and the rearrange_data comment still described ExaCA data as y -> x -> z
(fastest to slowest). ExaCA now writes x -> y -> z, which is what rearrange_data
already assumes, so only the documentation changes.
Remaining audit items that can't be fixed here: quick-xml advisories come in through
polars 0.55 -> object_store 0.13 (optional cloud storage, never compiled for this
crate, and voxel_coarsen does not commit a lockfile), and lexical unsoundness warnings
come from the rust_data_reader git dependency.
- xtal_light_up: built the old (0.21) and new (0.29) versions as Python modules and ran
all seven functions on the same random inputs. All 11 output arrays are bit-identical,
and the non-contiguous input error path still raises a Python exception.
- voxel_coarsen: default, python and polars features all build. On synthetic ExaCA-style
files (x -> y -> z order, including a box that is not square in x-y), the old build,
the new default build, the new polars build and the Python module all return the same
box size and coarsened data, and each block gets the expected majority grain.
Coarsening by a size that doesn't divide the box raises a Python exception as before.
… in-situ light-up
The xtal_light_up post-processing module (Rust with PyO3 bindings) could only handle
cubic crystals: the in-fiber check used the 24 cubic symmetry operators and a
reciprocal lattice built from a single lattice parameter, and the Taylor factor
rebuilt the effective plastic deformation rate from 12 hard-coded FCC slip systems.
ExaConstit's in-situ light-up (LightUp / LatticeTypeGeneral in mechanics_lightup.cpp)
supports all eight crystal systems, so the module now mirrors it and gives the same
results for every one of them.
New lattice.rs mirrors GetSymmetryGroups and LatticeTypeGeneral: the same eight
systems and names as the laue_type option (cubic, hexagonal, trigonal, rhombohedral,
tetragonal, orthorhombic, monoclinic, triclinic), the same symmetry operators, and the
same lattice parameter mapping and validation as lattice_params (parameter count per
system, positive values, angles in radians), along with a general reciprocal lattice
calculation.
calc_within_fibers gains an optional lattice_type = "cubic" argument, and
lattice_params now takes either a single value, as before, or a list. Existing calls
like calc_within_fibers(q, s, hkl, 3.60, tol, out) work unchanged. The wrong number of
parameters, non-positive parameters or an unknown lattice type raise an error.
The Taylor factor is now sum(|shear rates|) / equivalent plastic strain rate, using
ExaCMech's own equivalent plastic strain rate ("Equivalent Plastic Strain Rate"
output) the same way the in-situ light-up does, with the same 1e-14 cutoff. It works
for any number of slip systems. The FCC slip system / Schmid tensor code used to
rebuild that rate is removed rather than kept as a fallback, since it duplicated what
ExaCMech already outputs.
- calc_taylor_factors(shear_rates, eq_pl_strain_rate, taylor_factors)
- calc_taylor_factors_lattice_fiber(shear_rates, eq_pl_strain_rate, taylor_factors,
lattice_eq_pl_strain_rate, elem_vols, in_fibers, per_rank_update)
Calls with the old argument order fail with an error instead of giving wrong numbers,
and mismatched array sizes are checked.
- light_up_py/fiber_calcs_rank.py reads "Equivalent Plastic Strain Rate", uses the new
Taylor factor calls, and has new -lattice_type, -lattice_params and
-distance_tolerance_deg options. The defaults (cubic, 3.60, 5 degrees) reproduce
the previous behavior.
- pyxtallu.py drops its hand-built FCC Schmid tensor / equivalent plastic strain rate
calculation and passes an equivalent plastic strain rate in directly.
The trigonal list in GetSymmetryGroups repeated the a2 twofold axis where the a3 one
should be ({PI, -0.5, SQRT3_2, 0.0} twice instead of -SQRT3_2 for a3, as in the HEXRD
list it was taken from). That left 5 distinct operators that don't form a closed
group, so one symmetrically equivalent direction was never checked. This only changes
in-situ light-up results for trigonal materials.
- Doc comments throughout the Rust code; those on the Python-facing functions (array
shapes, component order, units, how outputs are filled, per_rank_update) also show
up in Python through help().
- Docstrings for both Python scripts, setup.py and the package __init__.py.
fiber_calcs_rank.py documents the ExaConstit settings it needs.
- New README.md covering building / installing, which ExaConstit outputs are needed
(adios2 visualization output, the six projected fields and the MPI rank count), the
crystal system table, an API summary, an example and the development workflow.
- src/options.toml listed the volume / centroid projections as "Volume" and
"Centroid"; the actual output names are "Element Volumes" and "Element Centroids".
- cargo clippy is clean (with -D warnings) after fixing the 26 existing warnings
(unused lifetimes, chunks_exact with constant sizes -> as_chunks, index loops, a
legacy numeric constant, doc comment placement), and the crate is run through
cargo fmt, which accounts for most of the diff in mod.rs. as_chunks needs Rust 1.88,
now recorded as rust-version in Cargo.toml (PyO3 already needed 1.83).
- New Rust unit tests: every symmetry group is closed under composition with no
duplicate operators, the reciprocal lattice is the exact dual of the direct lattice
for all eight systems, cubic B is I / a, and parameter validation.
- Cross-checked against ExaConstit itself: a small C++ program linked against the
ExaConstit library, using LatticeTypeGeneral and the same in-fiber loop as
LightUp::CalculateInFibers, compared with the Rust module over 3000 random
orientations x 8 HKLs for each crystal system. 0 mismatches for all eight.
- Cubic results are bit-identical to the previous version (lattice strains,
directional stiffness, in-fiber flags), with lattice_params given as a single value
or a list.
- Taylor factors match an independent NumPy calculation for 12, 18 and 48 slip
systems and both per_rank_update settings; the error cases raise.
- fiber_calcs_rank.py run end to end against a mock ADIOS2 reader (hexagonal, 18 slip
systems), with its saved fiber averages matching NumPy. Not yet run against a real
ADIOS2 .bp output.
- ExaConstit rebuilt with the trigonal fix.
Member
Author
|
Went through a bit of a bug-fixing spree and also fixing a few other updates on the Rust side of things to make sure some of the post-processing scripts had up to date packages and weren't hitting some warnings from |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Port of SNLS's trust region dogleg solver over to ExaConstit - I had Claude largely do this as this would have been a rather simple port as the math from both codes maps fairly easily, and I needed to test a wide range of new solvers to try and make some head ways for some tougher material modelling cases... Claude honestly did a good job and even did a complete port of the BBar method to the PA, so we fully support that code over there now as well...
Also I'm including a number of fun bug fixes that I've accumulated since the last release. Some of these were quite non-trivial and I spent quite a bit of time trying to nail down what was causing things to awol in some problems. These were largely MPI stalls or new GPU bugs found when moving to new versions of MFEM. For these, I had thrown codex at them after figuring out largely where the bugs might be but not be able to work out the exact locations when they were happening at scale. It was able to automate a lot of the MPI logging aspect and work out where things were stalling or diverging much faster than I would have been able to type. So definitely a nice debugger tool to use going forward...