Consistent tangent

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## latest
 
+- Add consistent homogenized tangent from the materials' own tangents (NEML2 and native models) [#163](https://github.com/DataAnalyticsEngineering/FANS/pull/163)
 - 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)
 

CMakeLists.txt

@@ -156,7 +156,8 @@ target_include_directories(FANS_main PRIVATE "${PROJECT_BINARY_DIR}/include")
 # SOURCES
 # ##############################################################################
 
-target_sources(FANS_FANS PRIVATE src/reader.cpp src/plugin_material.cpp)
+target_sources(FANS_FANS PRIVATE src/reader.cpp src/plugin_material.cpp
+                                 src/tangent.cpp)
 target_sources(FANS_main PRIVATE src/main.cpp)
 
 # ##############################################################################

docs/NEML2.md

@@ -27,10 +27,10 @@ This gives `test/FANS_neml2`. Without pixi, configure with `-DFANS_NEML2=ON` in
 ## Compile a model
 
 ```bash
-pixi run -e dev-neml2 neml2-compile model.i "cpu cuda"
+pixi run -e dev-neml2 neml2-compile model.i --model model --device cpu cuda --dtype float64 --output-dir compiled_models -d state/S:forces/E
 ```
 
-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"`.
+This compiles the block `[model]` of `model.i` for both devices into `compiled_models/model`. `--load a.py` (repeatable) imports a Python file the model needs. `-d flux:gradient` also compiles the model's tangent, which the `homogenized_tangent` result needs.
 
 ## Input file
 

include/MaterialManager.h

@@ -138,6 +138,13 @@ class MaterialManager {
         }
     }
 
+    // One linear model for all phases, e.g. the frozen tangent of src/tangent.cpp
+    MaterialManager(Matmodel<howmany, n_str> *model, int n_phases)
+        : phase_to_info(new MaterialInfo<howmany, n_str>[n_phases]), n_phases(n_phases), models{model}, all_stiffness(false)
+    {
+        std::fill_n(phase_to_info, n_phases, MaterialInfo<howmany, n_str>{model, nullptr, 0, true});
+    }
+
     ~MaterialManager()
     {
         for (auto *m : models)
@@ -221,10 +228,6 @@ class MaterialManager {
     {
         return phase_to_info[phase_id];
     }
-    inline size_t get_num_models() const
-    {
-        return models.size();
-    } // Get number of unique material models
     inline int get_num_phases() const
     {
         return n_phases;

include/fans_plugin.h

@@ -53,10 +53,11 @@ FANSPluginModel *fans_plugin_load(const char *config, char *msg, size_t msglen);
 int fans_plugin_set_table(FANSPluginModel *model, size_t field, int per_phase, size_t n_rows, const double *table,
                           char *err, size_t errlen);
 
-/* The step runs from time t_old to t. */
+/* The step runs from time t_old to t. Unless NULL, `tangent` gets the consistent
+   tangent d flux / d gradient, [n_points][n_str][n_str]. */
 int fans_plugin_evaluate(FANSPluginModel *model, size_t n_points, double t_old, double t, const double *gradient,
                          const int *voxel, const int *phase, const double *history_old, double *flux,
-                         double *history_new, char *err, size_t errlen);
+                         double *history_new, double *tangent, char *err, size_t errlen);
 
 void fans_plugin_free(FANSPluginModel *model);
 

include/material_models/large_strain/CompressibleNeoHookean.h

@@ -47,6 +47,26 @@ class CompressibleNeoHookean : public LargeStrainMechModel {
         return lambda[mat_index] * logJ * C_inv + mu[mat_index] * (Matrix3d::Identity() - C_inv);
     }
 
+    bool has_tangent() const override
+    {
+        return true;
+    }
+
+    // A = dP/dF of P = mu F + (lambda log(J) - mu) G with G = F^{-T}:
+    //   A_iJkL = mu d_ik d_JL + lambda G_iJ G_kL + (mu - lambda log(J)) G_iL G_kJ
+    void get_tangent(int i_eps, int mat_index, ptrdiff_t element_idx, Tangent A) override
+    {
+        const Matrix3d F = extract_F(i_eps);
+        const Matrix3d G = F.inverse().transpose();
+        const double   c = mu[mat_index] - lambda[mat_index] * log(F.determinant());
+
+        for (int i = 0; i < 3; ++i)
+            for (int J = 0; J < 3; ++J)
+                for (int k = 0; k < 3; ++k)
+                    for (int L = 0; L < 3; ++L)
+                        A(3 * i + J, 3 * k + L) = mu[mat_index] * (i == k) * (J == L) + lambda[mat_index] * G(i, J) * G(k, L) + c * G(i, L) * G(k, J);
+    }
+
     // At F = I the tangent is that of linear isotropic elasticity
     Matrix<double, 9, 9> get_reference_stiffness() override
     {

include/material_models/thermal/GBDiffusion.h

@@ -74,8 +74,8 @@ class GBDiffusion : public ThermalModel, public LinearModel<1, 3> {
             const int num_GB       = static_cast<int>(boundary_count);
             n_mat                  = num_crystals + num_GB;
 
-            const vector<double> grain_rot_matrices = Reader::ReadData<double>(reader.ms_filename, reader.MSGroup() + "rotation_matrices");
-            const vector<double> GB_normals         = Reader::ReadData<double>(reader.ms_filename, reader.MSGroup() + "GB_normals");
+            const vector<double> grain_rot_matrices = Reader::ReadData<double>(reader.ms_filename, reader.ms_group + "rotation_matrices");
+            const vector<double> GB_normals         = Reader::ReadData<double>(reader.ms_filename, reader.ms_group + "GB_normals");
             if (grain_rot_matrices.size() < 9 * size_t(num_crystals) || GB_normals.size() < 3 * size_t(n_mat))
                 throw std::runtime_error("rotation_matrices or GB_normals has too few entries");
 

include/matmodel.h

@@ -12,13 +12,9 @@ class Matmodel {
   public:
     static constexpr int num_str = n_str; // length of strain and stress
 
-    int    verbosity; //!< output verbosity
-    int    n_mat;     //!< Number of Materials
-    int    n_gp;      //!< Number of Gauss points (computed from FE_type)
-    string FE_type;   //!< Finite element type: "HEX8", "HEX8R", "BBAR"
-
-    double *strain; //!< Gradient
-    double *stress; //!< Flux
+    int    n_mat;   //!< Number of Materials
+    int    n_gp;    //!< Number of Gauss points (computed from FE_type)
+    string FE_type; //!< Finite element type: "HEX8", "HEX8R", "BBAR"
 
     Matmodel(const Reader &reader);
 
@@ -27,6 +23,22 @@ class Matmodel {
     void                                     getStrainStress(double *strain, double *stress, Matrix<double, howmany * 8, 1> &ue, int mat_index, ptrdiff_t element_idx);
     void                                     setGradient(vector<double> _g0);
 
+    // As getStrainStress, for the tangents C = d sigma / d eps of the converged step at the
+    // element's Gauss points, [n_gp][n_str][n_str]
+    using Tangent = Map<Matrix<double, n_str, n_str, RowMajor>>;
+    void getTangent(double *tangent, Matrix<double, howmany * 8, 1> &ue, int mat_index, ptrdiff_t element_idx)
+    {
+        eps.noalias() = B * ue + g0;
+        for (int i = 0; i < n_gp; ++i)
+            get_tangent(n_str * i, mat_index, element_idx, Tangent(tangent + i * n_str * n_str));
+    }
+    // A nonlinear model with a tangent: a batched one gives it in evaluate_batch,
+    // any other writes a get_tangent and says so here
+    virtual bool has_tangent() const
+    {
+        return wants_batch();
+    }
+
     // Accessors for internal Gauss point data (populated after getStrainStress call)
     inline const double *get_eps_data() const
     {
@@ -36,10 +48,6 @@ class Matmodel {
     {
         return sigma.data();
     }
-    inline int get_n_gp() const
-    {
-        return n_gp;
-    }
 
     virtual void postprocess(Solver<howmany, n_str> &solver, Reader &reader, int load_idx, int time_idx) {};
 
@@ -52,8 +60,10 @@ class Matmodel {
     {
         return false;
     }
+    // Unless null, tangent_gp gets the tangent d flux / d gradient of the converged
+    // step, [element][n_gp][n_str][n_str] (src/tangent.cpp).
     virtual void evaluate_batch(const vector<ptrdiff_t> &elems, const unsigned short *phase,
-                                const double *ue, double *sig_gp) {}
+                                const double *ue, double *sig_gp, double *tangent_gp = nullptr) {}
 
     // A model without the element stiffnesses of a LinearModel can still say that its
     // flux is linear in the gradient: then the linear CG applies.
@@ -87,6 +97,16 @@ class Matmodel {
     void                                       Construct_B();
 
     virtual void get_sigma(int i, int mat_index, ptrdiff_t element_idx) = 0;
+    // As get_sigma, for the tangent C = d sigma / d eps at that Gauss point. This one is that
+    // of a linear model, its stresses at unit strains; a nonlinear model writes its own.
+    virtual void get_tangent(int i, int mat_index, ptrdiff_t element_idx, Tangent C)
+    {
+        for (int j = 0; j < n_str; ++j) {
+            eps.segment<n_str>(i) = Matrix<double, n_str, 1>::Unit(j);
+            get_sigma(i, mat_index, element_idx);
+            C.col(j) = sigma.segment<n_str>(i);
+        }
+    }
 };
 
 template <int howmany, int n_str>

include/reader.h

@@ -25,9 +25,9 @@ class Reader {
     ~Reader();
 
     // contents of input file:
-    char             ms_filename[4096]{};    // Name of Micro-structure hdf5 file
-    char             ms_datasetname[4096]{}; // Absolute path of Micro-structure in hdf5 file
-    char             results_prefix[4096]{};
+    string           ms_filename;    // Name of Micro-structure hdf5 file
+    string           ms_datasetname; // Absolute path of Micro-structure in hdf5 file
+    string           ms_group;       // Its group, where the datasets that go with it are: "/sphere/32x32x32/" for "/sphere/32x32x32/ms"
     int              n_mat;
     json             inputJson; // Complete input JSON (for MaterialManager)
     json             materialProperties;
@@ -45,9 +45,8 @@ class Reader {
     string           strain_type{"small"}; // "small" (default) or "large"
     string           FE_type;              // "HEX8" (default), "HEX8R", or "BBAR"
     vector<string>   resultsToWrite;
-    char             results_filename[4096]{}; // Output HDF5 filename
-    char             dataset_name[8192]{};     // Base path for results in HDF5 file
-    hid_t            results_file_id = -1;     // Open HDF5 file handle for results
+    string           dataset_name;         // Base path for results in HDF5 file
+    hid_t            results_file_id = -1; // Open HDF5 file handle for results
 
     // contents of microstructure file:
     vector<int>     dims;
@@ -66,12 +65,10 @@ class Reader {
     ptrdiff_t local_n1;
     ptrdiff_t local_1_start;
 
-    // void Setup(ptrdiff_t howmany);
     void ReadInputFile(const std::string &input_fn);
     void ReadMS(int hm);
     void FreeMS();
     void ComputeVolumeFractions();
-    // void ReadHDF5(char file_name[], char dset_name[]);
     void safe_create_group(hid_t file, const char *const name);
     void OpenResultsFile(const char *output_fn); // Open results file once
     void CloseResultsFile();                     // Explicitly close results file
@@ -94,10 +91,6 @@ class Reader {
     static std::vector<hsize_t> DataShape(const string &file, const string &dset_name);
     template <typename T>
     void ReadSlab(T *data, const std::vector<int> &extra_dims, const string &file, const string &dset_name) const;
-
-    // The microstructure's group, e.g. "/sphere/32x32x32/" for "/sphere/32x32x32/ms",
-    // where the datasets that go with it are
-    string MSGroup() const;
 };
 
 // Opens a dataset for reading; the caller closes it and file_id
@@ -262,9 +255,8 @@ void Reader::writeData(const char *fieldName, int load_idx, int time_idx, const
     if (std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) {
         return;
     }
-    char name[5096];
-    snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName);
-    WriteData(data, name, shape, rank);
+    const string name = dataset_name + "/load" + to_string(load_idx) + "/time_step" + to_string(time_idx) + "/" + fieldName;
+    WriteData(data, name.c_str(), shape, rank);
 }
 
 template <typename T>
@@ -273,9 +265,8 @@ void Reader::writeSlab(const char *fieldName, int load_idx, int time_idx, const
     if (std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) {
         return;
     }
-    char name[5096];
-    snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName);
-    WriteSlab(data, extra_dims, name);
+    const string name = dataset_name + "/load" + to_string(load_idx) + "/time_step" + to_string(time_idx) + "/" + fieldName;
+    WriteSlab(data, extra_dims, name.c_str());
 }
 
 #endif