Add NEML2 material plugin

.gitignore

@@ -212,7 +212,15 @@ test/input_files/**/*.json
 !test_FiniteStrainJ2Plasticity.json
 !test_MixedBCs_LargeStrain.json
 
+output/
 # pixi environments
 .pixi
 *.egg-info
 *.conda
+
+# NEML2 artifacts
+build-neml2/
+FANS_neml2
+test/NEML2/models/*/
+test/NEML2/models/*_aoti.i
+test/NEML2/

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## latest
 
+- Add NEML2 material support [#162](https://github.com/DataAnalyticsEngineering/FANS/pull/162)
 - Account for grain orientations in GBDiffusion material model [#161](https://github.com/DataAnalyticsEngineering/FANS/pull/161)
 
 ## v0.8.1

CMakeLists.txt

@@ -12,6 +12,7 @@ project(
 include(GNUInstallDirs)
 include(CMakePackageConfigHelpers)
 
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
@@ -127,8 +128,9 @@ add_executable(FANS_main)
 add_custom_command(
   TARGET FANS_main
   POST_BUILD
-  COMMAND ${CMAKE_COMMAND} -E create_symlink $<TARGET_FILE:FANS_main>
-          ${CMAKE_CURRENT_SOURCE_DIR}/test/$<TARGET_FILE_NAME:FANS_main>
+  COMMAND
+    ${CMAKE_COMMAND} -E create_symlink $<TARGET_FILE:FANS_main>
+    ${CMAKE_CURRENT_SOURCE_DIR}/test/$<IF:$<TARGET_EXISTS:fans_neml2>,FANS_neml2,$<TARGET_FILE_NAME:FANS_main>>
   COMMENT
     "Create a symlink for FANS executable to ${CMAKE_CURRENT_SOURCE_DIR}/test/")
 
@@ -154,14 +156,14 @@ target_include_directories(FANS_main PRIVATE "${PROJECT_BINARY_DIR}/include")
 # SOURCES
 # ##############################################################################
 
-target_sources(FANS_FANS PRIVATE src/reader.cpp)
+target_sources(FANS_FANS PRIVATE src/reader.cpp src/plugin_material.cpp)
 target_sources(FANS_main PRIVATE src/main.cpp)
 
 # ##############################################################################
 # LINKING
 # ##############################################################################
 
-target_link_libraries(FANS_FANS PRIVATE m)
+target_link_libraries(FANS_FANS PRIVATE m ${CMAKE_DL_LIBS})
 
 target_link_libraries(FANS_FANS PUBLIC HDF5::HDF5)
 if(HDF5_IS_PARALLEL)
@@ -230,3 +232,28 @@ if(FANS_ENABLE_TESTING)
   enable_testing()
   add_subdirectory(test)
 endif()
+
+# ##############################################################################
+# NEML2 MATERIAL PLUGIN
+# ##############################################################################
+
+option(FANS_NEML2 "Also build libfans_neml2.so, the NEML2 material plugin" OFF)
+if(FANS_NEML2)
+  find_package(Python REQUIRED COMPONENTS Interpreter)
+  execute_process(
+    COMMAND ${Python_EXECUTABLE} -c "import neml2; print(neml2.__path__[0])"
+    OUTPUT_VARIABLE neml2_dir
+    OUTPUT_STRIP_TRAILING_WHITESPACE COMMAND_ERROR_IS_FATAL ANY)
+  get_filename_component(torch_prefix "${Python_EXECUTABLE}/../.." ABSOLUTE)
+  set(ENV{torch_ROOT} ${torch_prefix})
+  find_package(neml2 CONFIG REQUIRED PATHS ${neml2_dir})
+
+  add_library(fans_neml2 SHARED src/neml2_plugin.cpp)
+  target_include_directories(fans_neml2 PRIVATE include)
+  target_link_libraries(fans_neml2 PRIVATE neml2::aoti
+                                           nlohmann_json::nlohmann_json)
+  set_target_properties(
+    fans_neml2 PROPERTIES CXX_STANDARD 20
+                          INSTALL_RPATH "${neml2_dir}/lib;${torch_prefix}/lib")
+  install(TARGETS fans_neml2 LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
+endif()

README.md

@@ -197,6 +197,8 @@ cd ../test
 | `FANS_BUILD_STATIC` | Build static library | `OFF` |
 | `CMAKE_INSTALL_PREFIX` | Installation directory | System default |
 | `FANS_ENABLE_SANITIZERS` | Enable runtime sanitizers (AddressSanitizer and LeakSanitizer) for memory debugging | `OFF` |
+| `FANS_ENABLE_NATIVE` | Optimise for the building machine's CPU (`-march=native`) | `OFF` |
+| `FANS_NEML2` | Build a FANS that runs [NEML2](https://github.com/applied-material-modeling/neml2) materials | `OFF` |
 
 ---
 
@@ -284,6 +286,8 @@ FANS requires a JSON input file specifying the problem parameters. Example input
     - `CompressibleNeoHookean` for the compressible Neo-Hookean material model.
     - `FiniteStrainJ2Plasticity` for rate-independent finite-strain J2 plasticity with linear isotropic hardening.
 
+    - `NEML2` for a thermal, small- or large-strain material compiled from a [NEML2](https://github.com/applied-material-modeling/neml2) model with `neml2-compile` (needs `FANS_NEML2`), see [docs/NEML2.md](docs/NEML2.md).
+
   - `material_properties`: Material parameters specific to the chosen model. Properties are defined as arrays, where each element corresponds to one of the phases listed in the `phases` array.
 
 ### Solver settings

docs/NEML2.md

@@ -0,0 +1,60 @@
+# NEML2 materials in FANS
+
+Any [NEML2](https://github.com/applied-material-modeling/neml2) material model can be the material of a phase: `"matmodel": "NEML2"`.
+
+## Capabilities
+
+- Thermal, small-strain and large-strain problems: the model maps the temperature gradient to the heat flux, the strain to the stress, or the deformation gradient to the first Piola stress.
+- Any model `neml2-compile` can compile, e.g. from NEML2's own building blocks (J2 or crystal plasticity with an implicit update).
+- Your own model is a few lines of [PyTorch](https://pytorch.org/), automatic differentiation included (e.g. a hyperelastic law from just its strain energy): once compiled, it runs in FANS as is.
+- CPU or GPU. On a GPU all MPI ranks share the card.
+- History variables are found, stored and written by FANS (`"internal_variables"` in `results`).
+- Further model inputs ("fields", e.g. a grain orientation) come from the microstructure file, per voxel or per phase.
+- Time-dependent models get the time step of the input file (`time_step`).
+- Native and NEML2 materials can be mixed in one simulation.
+- NEML2 and libtorch never reach FANS itself: the material is a plugin library, `libfans_neml2.so`, loaded at run time.
+
+On the CPU, expect a NEML2 material to run about 2 to 6 times slower than the same model written natively in FANS.
+
+## Build
+
+```bash
+pixi run -e dev-neml2 build-fans_neml2
+```
+
+This gives `test/FANS_neml2`. Without pixi, configure with `-DFANS_NEML2=ON` in a Python environment that has `neml2` (3.1 or newer) and PyTorch installed.
+
+## Compile a model
+
+```bash
+pixi run -e dev-neml2 neml2-compile model.i "cpu cuda"
+```
+
+This compiles the model block named like the file (`[model]`) for both devices into `compiled_models/model`, next to `model.i`. A third argument lists Python files the model needs, e.g. `"a.py b.py"`.
+
+## Input file
+
+```json
+{
+    "phases": [0, 1],
+    "matmodel": "NEML2",
+    "material_properties": {
+        "artifact": "compiled_models/model",
+        "gradient": "forces/E",
+        "flux": "state/S",
+        "device": "cuda",
+        "fields": {"orientation": "rotation_matrices"}
+    }
+}
+```
+
+| Property | Meaning | Default |
+| --- | --- | --- |
+| `artifact` | Folder of the compiled model | required |
+| `gradient`, `flux` | Names of the model's gradient input and flux output | required |
+| `device` | `"cpu"`, `"cuda"`, ... | `"cpu"` |
+| `batch_size` | Material points per evaluation | 1024 on the CPU, 65536 on a GPU |
+| `fields` | For every further model input, the dataset to read it from: next to the microstructure, or an absolute path; shaped `[Z][Y][X][...]` (per voxel) or `[n_phase][...]` (per phase) | required if the model has such inputs |
+| `linear` | `true` if the flux is linear in the gradient: FANS then solves with its linear CG. | `false` |
+
+The input file also needs a top-level `reference_material`, the reference stiffness of the fundamental solution.

include/LargeStrainMechModel.h

@@ -84,100 +84,21 @@ class LargeStrainMechModel : public Matmodel<3, 9> {
     {
         throw std::logic_error("Large-strain material model must override compute_S() or get_sigma().");
     }
-    virtual Matrix<double, 6, 6> compute_material_tangent(const Matrix3d &F, int mat_index, ptrdiff_t element_idx, int i)
-    {
-        throw std::logic_error("Algorithmic material tangent is not implemented for this large-strain model.");
-    }
 
     /**
-     * @brief Convert material tangent C (dS/dE) to spatial tangent A (dP/dF)
-     *
-     * Computes the spatial tangent dP/dF from the material tangent dS/dE.
-     *
-     * The full expression is:
-     *   dP_iJ/dF_kL = δ_ik S_LJ + F_iM * dS_MJ/dE_PQ * dE_PQ/dF_kL
-     * where:
-     *   dE_PQ/dF_kL = 0.5 * (F_kP δ_QL + F_kQ δ_PL)
-     *
-     * @param F Deformation gradient (3×3)
-     * @param S 2nd Piola-Kirchhoff stress (3×3 symmetric)
-     * @param C_mandel Material tangent dS/dE in Mandel notation (6×6)
-     * @return Spatial tangent A = dP/dF in full notation (9×9)
+     * @brief Reference stiffness dP/dF at F = I of an isotropic material with Lame
+     * constants lambda and mu, as a 9x9 matrix on row-major F and P:
+     *   A_iJkL = lambda d_iJ d_kL + mu (d_ik d_JL + d_iL d_Jk)
      */
-    inline Matrix<double, 9, 9> compute_spatial_tangent(const Matrix3d             &F,
-                                                        const Matrix3d             &S,
-                                                        const Matrix<double, 6, 6> &C_mandel) const
+    static Matrix<double, 9, 9> isotropic_reference_stiffness(double lambda, double mu)
     {
         Matrix<double, 9, 9> A = Matrix<double, 9, 9>::Zero();
-
-        // Mandel to tensor index mapping
-        int mandel_to_ij[6][2] = {{0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}};
-
-        // Compute A_iJkL = dP_iJ/dF_kL
-        for (int i = 0; i < 3; ++i) {
-            for (int J = 0; J < 3; ++J) {
-                int row = 3 * i + J; // Index for P_iJ
-
-                for (int k = 0; k < 3; ++k) {
-                    for (int L = 0; L < 3; ++L) {
-                        int col = 3 * k + L; // Index for F_kL
-
-                        // First term: δ_ik S_LJ
-                        if (i == k) {
-                            A(row, col) += S(L, J);
-                        }
-
-                        // Second term: F_iM * C_MJPQ * dE_PQ/dF_kL
-                        for (int M = 0; M < 3; ++M) {
-                            // Find MJ in Mandel notation
-                            int MJ_mandel = -1;
-                            for (int idx = 0; idx < 6; ++idx) {
-                                if ((mandel_to_ij[idx][0] == M && mandel_to_ij[idx][1] == J) ||
-                                    (mandel_to_ij[idx][0] == J && mandel_to_ij[idx][1] == M)) {
-                                    MJ_mandel = idx;
-                                    break;
-                                }
-                            }
-                            if (MJ_mandel < 0)
-                                continue;
-
-                            for (int P = 0; P < 3; ++P) {
-                                for (int Q = P; Q < 3; ++Q) { // Q >= P due to symmetry
-                                    // Find PQ in Mandel notation
-                                    int PQ_mandel = -1;
-                                    for (int idx = 0; idx < 6; ++idx) {
-                                        if ((mandel_to_ij[idx][0] == P && mandel_to_ij[idx][1] == Q) ||
-                                            (mandel_to_ij[idx][0] == Q && mandel_to_ij[idx][1] == P)) {
-                                            PQ_mandel = idx;
-                                            break;
-                                        }
-                                    }
-                                    if (PQ_mandel < 0)
-                                        continue;
-
-                                    // Get C_MJPQ and account for Mandel factors
-                                    double C_val = C_mandel(MJ_mandel, PQ_mandel);
-                                    if (MJ_mandel >= 3)
-                                        C_val /= sqrt(2.0);
-                                    if (PQ_mandel >= 3)
-                                        C_val /= sqrt(2.0);
-
-                                    // Compute dE_PQ/dF_kL = 0.5 * (F_kP δ_QL + F_kQ δ_PL)
-                                    double dE_dF = 0.0;
-                                    if (Q == L)
-                                        dE_dF += 0.5 * F(k, P);
-                                    if (P == L)
-                                        dE_dF += 0.5 * F(k, Q);
-
-                                    A(row, col) += F(i, M) * C_val * dE_dF;
-                                }
-                            }
-                        }
-                    }
-                }
+        for (int i = 0; i < 3; ++i)
+            for (int j = 0; j < 3; ++j) {
+                A(3 * i + i, 3 * j + j) += lambda; // lambda d_iJ d_kL
+                A(3 * i + j, 3 * i + j) += mu;     // mu d_ik d_JL
+                A(3 * i + j, 3 * j + i) += mu;     // mu d_iL d_Jk
             }
-        }
-
         return A;
     }