Nonlinear mechanics with SUNDIALS#469
Conversation
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request integrates the SUNDIALS KINSOL nonlinear solver into the mechanical physics module and updates the build environment across Docker, Nix, and CMake. The review identifies a package naming error in the Dockerfile and several implementation gaps in the solver logic, specifically the need to delegate residual calculations to the operator for true nonlinear support, implement the Jacobian setup callback, and verify the solver's convergence status.
| python3-dev \ | ||
| git \ | ||
| lcov \ | ||
| sundials \ |
There was a problem hiding this comment.
The package name for SUNDIALS development headers in Ubuntu is libsundials-dev. Using sundials will likely result in an 'Unable to locate package' error, causing the Docker build to fail.
libsundials-dev \
| nonlinear_solver.residual = [&](const TrilinosVectorType &evaluation_point, | ||
| TrilinosVectorType &residual) | ||
| { | ||
| _mechanical_operator->system_matrix().vmult(residual, evaluation_point); | ||
| residual -= rhs_device; | ||
| }; |
There was a problem hiding this comment.
The residual calculation is currently hardcoded to use the linear system matrix: _mechanical_operator->system_matrix().vmult(...). This effectively restricts the nonlinear solver to solving a linear system, which adds unnecessary overhead. To properly support nonlinear mechanics as intended, the residual calculation should be delegated to the MechanicalOperator (e.g., via a residual() method) to allow for the inclusion of nonlinear terms.
| nonlinear_solver.setup_jacobian = | ||
| [&](const TrilinosVectorType & /*current_u*/, | ||
| const TrilinosVectorType & /*current_f*/) {}; |
There was a problem hiding this comment.
The setup_jacobian callback is currently empty. For nonlinear problems, KINSOL relies on this callback to update the Jacobian (tangent stiffness matrix) when necessary. Without an implementation that updates the matrix based on the current solution state, the solver is limited to a Modified Newton method using the initial linear stiffness, which may fail to converge or be very inefficient for significant nonlinearities.
| _mechanical_operator->preconditioner()); | ||
| }; | ||
|
|
||
| nonlinear_solver.solve(displacement); |
There was a problem hiding this comment.
The return code of nonlinear_solver.solve(displacement) is ignored. If KINSOL fails to converge, the simulation will proceed with an incorrect solution. It is important to check the return status and handle failures (e.g., by throwing an exception or logging a critical error).
|
With the linear operator currently used, I'm seeing for the diff --git a/source/MechanicalPhysics.cc b/source/MechanicalPhysics.cc
index f57fb4c..7c7027d 100644
--- a/source/MechanicalPhysics.cc
+++ b/source/MechanicalPhysics.cc
@@ -510,9 +510,11 @@ MechanicalPhysics<dim, n_materials, p_order, MaterialStates,
dealii::SolverCG<TrilinosVectorType> cg(solver_control);
cg.solve(_mechanical_operator->system_matrix(), dst, rhs,
_mechanical_operator->preconditioner());
+ std::cout << "linear_iterations: " << solver_control.last_step() << std::endl;
};
- nonlinear_solver.solve(displacement);
+ auto nonlinear_iterations = nonlinear_solver.solve(displacement);
+ std::cout << "nonlinear_iterations: " << nonlinear_iterations << std::endl;
rw_vector.import_elements(displacement, dealii::VectorOperation::insert);
dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>the following |
|
I need to rebuild the base image for this PR to run. Modifying the Dockerfile does not trigger a rebuild of the image |
Right but I confirmed that the updated Dockerfile works. |
| wget --quiet ${P4EST_URL} --output-document=${P4EST_ARCHIVE} && \ | ||
| mkdir -p ${P4EST_SOURCE_DIR} && \ | ||
| cd ${P4EST_SOURCE_DIR} && \ | ||
| wget --quiet https://raw.githubusercontent.com/dealii/dealii/master/doc/external-libs/p4est-setup.sh && \ |
There was a problem hiding this comment.
I had trouble accessing that URL from ORNL networks.
Shouldn't this converge in one iteration since the problem is linear? |
KINSOL chooses the step size in the update as well as the tolerance for every iteration step. If the step size was 1 and the solver tolerance was the same as the final tolerance, we should converge in one step. |
|
Does this PR just need a new Docker image and a fix for the Nix build? |
Yes. Of course, solving the Mechanical problem also is (unnecessarily) more expensive so the question is if we care about a shortcut if we know that the constitutive law is linear. |
|
We probably want a shortcut. My understanding is that the people at MDF first use a linear model to narrow down the build parameters and then, they use a more expensive model to refine the parameters. So we want the linear to build as fast as possible since that's what they run the most. |
|
I created new two images |
masterleinad
force-pushed
the
nonlinear_mechanics_sundials
branch
from
99e47fc to
769e3cb
Compare
| python3-dev \ | ||
| git \ | ||
| lcov \ | ||
| libsundials-dev \ |
There was a problem hiding this comment.
This pulls MPI, so we cannot install our own MPI anymore.
There was a problem hiding this comment.
OK. Let me try building SUNDIALS explicitly in the image then.
There was a problem hiding this comment.
We can also just the package manager as much as possible and stop installing our own library. I don't what's better.
One reason I had to install everything by hand, was to show people how to install the software stack on their own machine. As far as I know, nobody does this, instead people are using the Nix build.
ARG BASE=nvidia/cuda:12.6.2-devel-ubuntu24.04
FROM $BASE
ARG N_PROCS=8
ARG CUDA=ON
ENV DEBIAN_FRONTEND=noninteractive
RUN apt-get update && apt-get upgrade -y && apt-get install -y \
cmake \
gcc \
gfortran \
clang-format \
build-essential \
wget \
curl \
libboost-all-dev \
libcurl4-gnutls-dev \
bison \
python3-dev \
git \
lcov \
zlib1g-dev \
libopenblas-dev \
libdw-dev \
libpfm4-dev \
libunwind-dev \
libgmp-dev \
openmpi-bin \
libopenmpi-dev \
libp4est-dev \
libhdf5-dev \
libnetcdf-dev \
libsundials-dev \
&& \
apt-get clean && rm -rf /var/lib/apt/lists/*
ENV PREFIX=/scratch
ENV ARCHIVE_DIR=${PREFIX}/archive
ENV SOURCE_DIR=${PREFIX}/source
ENV BUILD_DIR=${PREFIX}/build
ENV INSTALL_DIR=/opt
RUN mkdir -p ${PREFIX} && \
cd ${PREFIX} && \
mkdir archive && \
mkdir source && \
mkdir build
ENV OMPI_ALLOW_RUN_AS_ROOT=1
ENV OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1
ENV OMPI_MCA_btl_vader_single_copy_mechanism=none
ENV OMPI_MCA_btl_base_warn_component_unused='0'
ENV OMPI_CC=cc
ENV OMPI_CXX=c++
# Install Trilinos 14.4.0
RUN export TRILINOS_HASH=975307431d60d0859ebaa27c9169cbb1d4287513 && \
export TRILINOS_URL=https://github.com/trilinos/Trilinos/archive/${TRILINOS_HASH}.tar.gz && \
export TRILINOS_ARCHIVE=${ARCHIVE_DIR}/trilinos.tar.gz && \
export TRILINOS_SOURCE_DIR=${SOURCE_DIR}/trilinos && \
export TRILINOS_BUILD_DIR=${BUILD_DIR}/trilinos && \
export TRILINOS_INSTALL_DIR=${INSTALL_DIR}/trilinos && \
wget ${TRILINOS_URL} --output-document=${TRILINOS_ARCHIVE} && \
mkdir -p ${TRILINOS_SOURCE_DIR} && \
tar -xf ${TRILINOS_ARCHIVE} -C ${TRILINOS_SOURCE_DIR} --strip-components=1 && \
mkdir -p ${TRILINOS_INSTALL_DIR} && \
cp -r ${TRILINOS_SOURCE_DIR}/packages/kokkos/bin ${TRILINOS_INSTALL_DIR} && \
export OMPI_CXX=${TRILINOS_INSTALL_DIR}/bin/nvcc_wrapper && \
mkdir ${TRILINOS_BUILD_DIR} && \
cd ${TRILINOS_BUILD_DIR} && \
cmake -DCMAKE_BUILD_TYPE=RELEASE -DBUILD_SHARED_LIBS=ON \
-DTPL_ENABLE_MPI=ON \
-DTPL_ENABLE_BLAS=ON \
-DTPL_ENABLE_LAPACK=ON \
-DTPL_ENABLE_Boost=ON \
-DBoost_INCLUDE_DIRS=${BOOST_DIR}/include \
-DTPL_ENABLE_BoostLib=ON \
-DBoostLib_INCLUDE_DIRS=${BOOST_DIR}/include \
-DBoostLib_LIBRARY_DIRS=${BOOST_DIR}/lib \
-DTrilinos_ENABLE_ALL_PACKAGES=OFF \
-DTrilinos_ENABLE_ALL_OPTIONAL_PACKAGES=OFF \
-DTrilinos_ENABLE_TESTS=OFF \
-DTrilinos_ENABLE_EXAMPLES=OFF \
-DTrilinos_ENABLE_Amesos=ON \
-DTrilinos_ENABLE_AztecOO=ON \
-DTrilinos_ENABLE_Epetra=ON \
-DTrilinos_ENABLE_Ifpack=ON \
-DTrilinos_ENABLE_Kokkos=ON \
-DTrilinos_ENABLE_ML=ON \
-DTrilinos_ENABLE_SEACAS=ON \
-DTrilinos_ENABLE_EXPLICIT_INSTANTIATION=ON \
-DKokkos_ENABLE_CUDA=$CUDA \
-DKokkos_ENABLE_CUDA_CONSTEXPR=ON \
-DKokkos_ARCH_AMPERE80=ON \
-DBLAS_LIBRARY_NAMES="openblas" \
-DBLAS_LIBRARY_DIRS=${OPENBLAS_DIR}/lib \
-DLAPACK_LIBRARY_NAMES="openblas" \
-DLAPACK_LIBRARY_DIRS=${OPENBLAS_DIR}/lib \
-DTPL_ENABLE_Netcdf=ON \
-DNetcdf_INCLUDE_DIRS=$NETCDF_DIR/include \
-DNetcdf_LIBRARY_DIRS=$NETCDF_DIR/lib \
-DTPL_ENABLE_Matio=OFF \
-DTPL_ENABLE_X11=OFF \
-DCMAKE_INSTALL_PREFIX=${TRILINOS_INSTALL_DIR} \
${TRILINOS_SOURCE_DIR} && \
make -j${N_PROCS} install && \
rm -rf ${TRILINOS_ARCHIVE} && \
rm -rf ${TRILINOS_BUILD_DIR} && \
rm -rf ${TRILINOS_SOURCE_DIR}
ENV TRILINOS_DIR=${INSTALL_DIR}/trilinos
ENV OMPI_CXX=${TRILINOS_DIR}/bin/nvcc_wrapper
# Install Adiak
RUN export ADIAK_VERSION=0.4.0 && \
export ADIAK_URL=https://github.com/LLNL/Adiak/releases/download/v${ADIAK_VERSION}/adiak-${ADIAK_VERSION}.tar.gz && \
export ADIAK_ARCHIVE=${ARCHIVE_DIR}/adiak.tar.gz && \
export ADIAK_SOURCE_DIR=${SOURCE_DIR}/adiak && \
export ADIAK_BUILD_DIR=${BUILD_DIR}/adiak && \
export ADIAK_INSTALL_DIR=${INSTALL_DIR}/adiak && \
wget ${ADIAK_URL} --output-document=${ADIAK_ARCHIVE} && \
mkdir -p ${ADIAK_SOURCE_DIR} && \
tar -xf ${ADIAK_ARCHIVE} -C ${ADIAK_SOURCE_DIR} --strip-components=1 && \
mkdir -p ${ADIAK_BUILD_DIR} && \
cd ${ADIAK_BUILD_DIR} && \
cmake -DCMAKE_BUILD_TYPE=Release \
-DBUILD_SHARED_LIBS=ON \
-DCMAKE_INSTALL_PREFIX=${ADIAK_INSTALL_DIR} \
${ADIAK_SOURCE_DIR} && \
make -j${N_PROCS} install && \
rm -rf ${ADIAK_ARCHIVE} && \
rm -rf ${ADIAK_BUILD_DIR} && \
rm -rf ${ADIAK_SOURCE_DIR}
ENV ADIAK_DIR=${INSTALL_DIR}/adiak
# Install Caliper 2.11.0
RUN export CALIPER_VERSION=2.11.0 && \
export CALIPER_URL=https://github.com/LLNL/Caliper/archive/refs/tags/v${CALIPER_VERSION}.tar.gz && \
export CALIPER_ARCHIVE=${ARCHIVE_DIR}/caliper.tar.gz && \
export CALIPER_SOURCE_DIR=${SOURCE_DIR}/caliper && \
export CALIPER_BUILD_DIR=${BUILD_DIR}/caliper && \
export CALIPER_INSTALL_DIR=${INSTALL_DIR}/caliper && \
wget ${CALIPER_URL} --output-document=${CALIPER_ARCHIVE} && \
mkdir -p ${CALIPER_SOURCE_DIR} && \
tar -xf ${CALIPER_ARCHIVE} -C ${CALIPER_SOURCE_DIR} --strip-components=1 && \
mkdir -p ${CALIPER_BUILD_DIR} && \
cd ${CALIPER_BUILD_DIR} && \
cmake -DCMAKE_BUILD_TYPE=Release \
-DBUILD_SHARED_LIBS=ON \
-DWITH_ADIAK=ON \
-DWITH_LIBDW=ON \
-DWITH_LIBPFM=ON \
-DWITH_LIBUNWIND=ON \
-DWITH_MPI=ON \
-DWITH_SAMPLER=ON \
-DCMAKE_INSTALL_PREFIX=${CALIPER_INSTALL_DIR} \
${CALIPER_SOURCE_DIR} && \
make -j${N_PROCS} install && \
rm -rf ${CALIPER_ARCHIVE} && \
rm -rf ${CALIPER_BUILD_DIR} && \
rm -rf ${CALIPER_SOURCE_DIR}
ENV CALIPER_DIR=${INSTALL_DIR}/caliper
# Install ArborX 1.5
RUN export ARBORX_HASH=e026f8237dbd0aaf2ddf2e3bf0e492071db854a2 && \
export ARBORX_URL=https://github.com/arborx/ArborX/archive/${ARBORX_HASH}.tar.gz && \
export ARBORX_ARCHIVE=${ARCHIVE_DIR}/arborx.tar.gz && \
export ARBORX_SOURCE_DIR=${SOURCE_DIR}/arborx && \
export ARBORX_BUILD_DIR=${BUILD_DIR}/arborx && \
export ARBORX_INSTALL_DIR=${INSTALL_DIR}/arborx && \
wget ${ARBORX_URL} --output-document=${ARBORX_ARCHIVE} && \
mkdir -p ${ARBORX_SOURCE_DIR} && \
tar -xf ${ARBORX_ARCHIVE} -C ${ARBORX_SOURCE_DIR} --strip-components=1 && \
mkdir -p ${ARBORX_BUILD_DIR} && \
cd ${ARBORX_BUILD_DIR} && \
cmake -DCMAKE_BUILD_TYPE=Release \
-DBUILD_SHARED_LIBS=ON \
-DCMAKE_CXX_COMPILER=${TRILINOS_DIR}/bin/nvcc_wrapper \
-DCMAKE_CXX_EXTENSIONS=OFF \
-DARBORX_ENABLE_MPI=ON \
-DCMAKE_PREFIX_PATH=${TRILINOS_DIR} \
-DCMAKE_INSTALL_PREFIX=${ARBORX_INSTALL_DIR} \
${ARBORX_SOURCE_DIR} && \
make -j${N_PROCS} install && \
rm -rf ${ARBORX_ARCHIVE} && \
rm -rf ${ARBORX_BUILD_DIR} && \
rm -rf ${ARBORX_SOURCE_DIR}
ENV ARBORX_DIR=${INSTALL_DIR}/arborx
# Install deal.II 9.6.2
RUN export DEAL_II_VERSION=9.6.2 && \
export DEAL_II_URL=https://github.com/dealii/dealii/releases/download/v${DEAL_II_VERSION}/dealii-${DEAL_II_VERSION}.tar.gz && \
export DEAL_II_ARCHIVE=${ARCHIVE_DIR}/dealii.tar.gz && \
export DEAL_II_SOURCE_DIR=${SOURCE_DIR}/dealii && \
export DEAL_II_BUILD_DIR=${BUILD_DIR}/dealii && \
export DEAL_II_INSTALL_DIR=${INSTALL_DIR}/dealii && \
wget ${DEAL_II_URL} --output-document=${DEAL_II_ARCHIVE} && \
mkdir -p ${DEAL_II_SOURCE_DIR} && \
tar -xf ${DEAL_II_ARCHIVE} -C ${DEAL_II_SOURCE_DIR} --strip-components=1 && \
mkdir -p ${DEAL_II_BUILD_DIR} && cd ${DEAL_II_BUILD_DIR} && \
cmake -DCMAKE_BUILD_TYPE=DebugRelease \
-DCMAKE_CXX_COMPILER=${TRILINOS_DIR}/bin/nvcc_wrapper \
-DCMAKE_CXX_STANDARD=17 \
-DCMAKE_CXX_EXTENSIONS=OFF \
-DDEAL_II_EARLY_DEPRECATIONS=ON \
-DDEAL_II_WITH_TBB=OFF \
-DDEAL_II_WITH_64BIT_INDICES=ON \
-DDEAL_II_WITH_COMPLEX_VALUES=OFF \
-DDEAL_II_WITH_MPI=ON \
-DDEAL_II_WITH_SUNDIALS=ON \
-DSUNDIALS_DIR=${SUNDIALS_DIR} \
-DDEAL_II_WITH_P4EST=ON \
-DP4EST_DIR=${P4EST_DIR} \
-DDEAL_II_WITH_ARBORX=ON \
-DARBORX_DIR=${ARBORX_DIR} \
-DDEAL_II_WITH_TRILINOS=ON \
-DTRILINOS_DIR=${TRILINOS_DIR} \
-DDEAL_II_TRILINOS_WITH_SEACAS=ON \
-DDEAL_II_COMPONENT_EXAMPLES=OFF \
-DCMAKE_INSTALL_PREFIX=${DEAL_II_INSTALL_DIR} \
${DEAL_II_SOURCE_DIR} && \
make -j${N_PROCS} install && \
rm -rf ${DEAL_II_ARCHIVE} && \
rm -rf ${DEAL_II_BUILD_DIR}
# We keep the source file for debugging purpose
ENV DEAL_II_DIR=${INSTALL_DIR}/dealiiseems to compile fine. |
masterleinad
force-pushed
the
nonlinear_mechanics_sundials
branch
from
f5dcd9a to
bcde6f7
Compare
|
@Rombur Let me know if you want me to separate the changes to the Dockerfile from this pull request. |
|
Yes, can you split this PR in two: one with Docker, Nix, and the build system and another one with the rest. |
masterleinad
force-pushed
the
nonlinear_mechanics_sundials
branch
from
bcde6f7 to
05184e3
Compare
|
|
||
| nonlinear_solver.setup_jacobian = | ||
| [&](const TrilinosVectorType & /*current_u*/, | ||
| const TrilinosVectorType & /*current_f*/) {}; |
There was a problem hiding this comment.
Why is this empty? Shouldn't you be using finite difference to approximate the jacobian or does kinsol switch to a fixed point iteration if you don't provide a jacobian?
There was a problem hiding this comment.
No, JFNK isn't accessible through deal.II dircetly AFAICT. I left it empty for now since it's redundant if the operator is linear.
There was a problem hiding this comment.
Looking more at the code, I am really confused with the function solve_with_jacobian. My understanding is that this function wants you to solve J(u) dst = rhs. Since we will use a Krylov method, we just need the application of the jacobian. So we need the matrix-vector multiplication J(u) v. We can approximate this using [F(u + ε·v) − F(u) ] / ε, where F(x) = 0 is the problem we are trying to solve. So you need to create a new operator such that vmult does [F(u + ε·v) − F(u) ] / ε by calling the _mechanical_operator twice (though you can save F(u) to skip one operator evaluation) . Am I totally off?
There was a problem hiding this comment.
For a linear operator
[F(u + ε·v) − F(u) ] / ε = Fv
so we can just use the existing oeprator. As soon as we have a real nonlinear operator, we obviously have to do something else. In particular, setup_jacobian needs to do something. I'm not convinced that we should do the approximation ourselves but rather leave it to KINSOL if necessary (which needs to be implented in deal.II) or compute the jacobian ourselves. At this point, it's not clear to me what kind if nonlinear models we are interested in and if it's feasible to compute the jacobian manually.
There was a problem hiding this comment.
We are not going to compute the jacobian manually any time soon. We need to show that we can support nonlinear constitutive law and then later, we can think about making the code more efficient.
If KINSOL has functions to help with JFNK, we should use them. We can use KINSOL outside of deal.II if it's too much work to integrate that in deal.II so close to the release.
|
@AshGannon this PR is relevant to you |
2 participants
This pull requests allows for the mechanical operator to be nonlinear relying on SUNDIALS/KINSOL for the nonlinear iteration.