Integrate pyzag bindings for neml2

.gitattributes

@@ -0,0 +1,2 @@
+
+*.ipynb	diff=jupyternotebook

.github/workflows/python.yml

@@ -47,10 +47,10 @@ jobs:
       - uses: actions/setup-python@v5
         with:
           python-version: "3.9"
-      - name: Install Python dependencies
-        run: pip install -r requirements.txt
       - name: Install PyTorch
         run: pip install torch==2.5.1
+      - name: Install Python dependencies
+        run: pip install -r requirements.txt
       - run: pip install -v .
       - run: pytest --junitxml=pytest.xml python/tests
       - name: Publish Test Results

include/neml2/models/LinearInterpolation.h

@@ -95,17 +95,6 @@ class LinearInterpolation : public Interpolation<T>
    */
   template <typename T2>
   T2 mask(const T2 & in, const torch::Tensor & m) const;
-
-  /// Batch shape of the interpolant, excluding the last dimension which is the interpolation axis
-  const TensorShape _interp_batch_sizes;
-  /// Starting abscissa of each interval
-  const Scalar & _X0;
-  /// Ending abscissa of each interval
-  const Scalar & _X1;
-  /// Starting ordinate of each interval
-  const T & _Y0;
-  /// Slope of each interval
-  const T & _slope;
 };
 
 template <typename T>

include/neml2/models/solid_mechanics/KocksMeckingIntercept.h

@@ -0,0 +1,53 @@
+// Copyright 2024, UChicago Argonne, LLC
+// All Rights Reserved
+// Software Name: NEML2 -- the New Engineering material Model Library, version 2
+// By: Argonne National Laboratory
+// OPEN SOURCE LICENSE (MIT)
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+// THE SOFTWARE.
+
+#pragma once
+
+#include "neml2/models/NonlinearParameter.h"
+
+namespace neml2
+{
+/**
+ * @brief A scalar-valued parameter defined by (C-B) / A
+ */
+class KocksMeckingIntercept : public NonlinearParameter<Scalar>
+{
+public:
+  static OptionSet expected_options();
+
+  KocksMeckingIntercept(const OptionSet & options);
+
+protected:
+  void set_value(bool out, bool dout_din, bool d2out_din2) override;
+
+  /// KM A
+  const Scalar & _A;
+
+  /// KM B
+  const Scalar & _B;
+
+  /// KM C
+  const Scalar & _C;
+};
+}

include/neml2/tensors/LabeledMatrix.h

@@ -44,10 +44,6 @@ class LabeledMatrix : public LabeledTensor<LabeledMatrix, 2>
                                 const LabeledAxis & axis,
                                 const torch::TensorOptions & options = default_tensor_options());
 
-  /// Fill another matrix into this matrix.
-  /// The item set of the other matrix must be a subset of this matrix's item set.
-  void fill(const LabeledMatrix & other, bool recursive = true);
-
   /// Chain rule product of two derivatives
   LabeledMatrix chain(const LabeledMatrix & other) const;
 };

include/neml2/tensors/LabeledTensor.h

@@ -160,6 +160,15 @@ class LabeledTensor
   /// Get a specific labeled axis
   const LabeledAxis & axis(Size i = 0) const { return *_axes[i]; }
 
+  /**
+   * @brief Fill with another LabeledTensor that matches this one on all but one axis
+   *
+   * @param Derived The LabeledTensor to fill from
+   * @param Size The common axis, default 0
+   * @param bool If true, fill recursively down subaxes
+   */
+  void fill(const Derived & other, Size odim = 0, bool recursive = true);
+
 protected:
   /// The tensor
   Tensor _tensor;

include/neml2/tensors/LabeledVector.h

@@ -41,10 +41,6 @@ class LabeledVector : public LabeledTensor<LabeledVector, 1>
 
   /// Slice the logically 1D tensor by a single sub-axis
   LabeledVector slice(const std::string & name) const;
-
-  /// Fill (override) another vector into this vector.
-  /// The item set of the other vector must be a subset of this vector's item set.
-  void fill(const LabeledVector & other, bool recursive = true);
 };
 
 namespace utils

python/CMakeLists.txt

@@ -43,6 +43,11 @@ install(FILES
   DESTINATION .
 )
 
+# pyzag subpackage
+configure_file(neml2/pyzag/__init__.py ${NEML2_BINARY_DIR}/python/neml2/pyzag/__init__.py COPYONLY)
+configure_file(neml2/pyzag/interface.py ${NEML2_BINARY_DIR}/python/neml2/pyzag/interface.py COPYONLY)
+install(DIRECTORY neml2/pyzag DESTINATION .)
+
 # Tests
 if(NEML2_TESTS)
   add_subdirectory(tests)

python/examples/README.md

@@ -0,0 +1,21 @@
+# Jupyter notebook examples of the NEML2/pyzag interface
+
+## Setup
+
+To run these notebooks install the packages given in `python/examples/requirements.txt`
+
+```
+pip install -r python/examples/requirements.txt
+```
+
+after building the NEML2 python package and installing the main set of python requirements in `requirements.txt` in the root directory.
+
+## Version control for Jupyter notebook examples
+
+We track these notebooks with the `nbdime` python tool, installed as part of the package requirements.  The first time you install this for use in a new repository you need to run
+
+```
+git-nbdiffdriver config --enable
+```
+
+to enable the hooks.

python/examples/demo_model.i

@@ -0,0 +1,210 @@
+[Tensors]
+  [T_controls]
+    type = Scalar
+    values = '300 347.36842105 394.73684211 442.10526316 489.47368421 536.84210526 584.21052632 631.57894737 678.94736842 726.31578947 773.68421053 821.05263158 868.42105263 915.78947368 963.15789474 1010.52631579 1057.89473684 1105.26315789 1152.63157895 1200'
+    batch_shape = '(20)'
+  []
+  [mu_values]
+    type = Scalar
+    values = '76670.48346056 75465.18012589 74314.80514263 73374.72880675 72651.54680595 71928.36480514 71120.75130575 70035.97830454 68951.20530333 67842.26597027 66399.97991161 65315.20691041 63884.85335476 62763.98151868 61373.80474086 59927.44073925 58481.07673765 56544.43551627 54599.93973483 52791.98473282'
+    batch_shape = '(20)'
+  []
+  [T_train]
+    type = Scalar
+    values = '300.0 600.0 900.0 1200.0'
+    batch_shape = '(4)'
+  []
+  [R_values]
+    type = Scalar
+    values = '300.0 200.0 100.0 50.0'
+    batch_shape = '(4)'
+  []
+  [d_values]
+    type = Scalar
+    values = '30.0 20.0 15.0 12.0'
+    batch_shape = '(4)'
+  []
+[]
+
+[Models]
+  [A]
+    type = ScalarConstantParameter
+    value = -8.679
+  []
+  [B]
+    type = ScalarConstantParameter
+    value = -0.744
+  []
+  [C]
+    type = ScalarConstantParameter
+    value = -5.41
+  []
+  [g0]
+    type = KocksMeckingIntercept
+    A = 'A'
+    B = 'B'
+    C = 'C'
+  []
+  [mandel_stress]
+    type = IsotropicMandelStress
+  []
+  [vonmises]
+    type = SR2Invariant
+    invariant_type = 'VONMISES'
+    tensor = 'state/internal/M'
+    invariant = 'state/internal/s'
+  []
+  [isoharden]
+    type = VoceIsotropicHardening
+    saturated_hardening = 'R'
+    saturation_rate = 'd'
+  []
+  [mu]
+    type = ScalarLinearInterpolation
+    argument = 'forces/T'
+    abscissa = 'T_controls'
+    ordinate = 'mu_values'
+  []
+  [R]
+    type = ScalarLinearInterpolation
+    argument = 'forces/T'
+    abscissa = 'T_train'
+    ordinate = 'R_values'
+  []
+  [d]
+    type = ScalarLinearInterpolation
+    argument = 'forces/T'
+    abscissa = 'T_train'
+    ordinate = 'd_values'
+  []
+  [ys]
+    type = KocksMeckingYieldStress
+    shear_modulus = 'mu'
+    C = 'C'
+  []
+  [yield]
+    type = YieldFunction
+    yield_stress = 'ys'
+    isotropic_hardening = 'state/internal/k'
+  []
+  [yield_zero]
+    type = YieldFunction
+    yield_stress = 0
+    isotropic_hardening = 'state/internal/k'
+    yield_function = 'state/internal/fp_alt'
+  []
+  [flow]
+    type = ComposedModel
+    models = 'vonmises yield'
+    automatic_nonlinear_parameter = false
+  []
+  [normality]
+    type = Normality
+    model = 'flow'
+    function = 'state/internal/fp'
+    from = 'state/internal/M state/internal/k'
+    to = 'state/internal/NM state/internal/Nk'
+  []
+  [ri_flowrate]
+    type = RateIndependentPlasticFlowConstraint
+    flow_rate = 'state/internal/gamma_rate_ri'
+  []
+  [km_sensitivity]
+    type = KocksMeckingRateSensitivity
+    A = 'A'
+    shear_modulus = 'mu'
+    k = 1.38064e-20
+    b = 2.474e-7
+  []
+  [km_viscosity]
+    type = KocksMeckingFlowViscosity
+    A = 'A'
+    B = 'B'
+    shear_modulus = 'mu'
+    k = 1.38064e-20
+    b = 2.474e-7
+    eps0 = 1e10
+  []
+  [rd_flowrate]
+    type = PerzynaPlasticFlowRate
+    reference_stress = 'km_viscosity'
+    exponent = 'km_sensitivity'
+    yield_function = 'state/internal/fp_alt'
+    flow_rate = 'state/internal/gamma_rate_rd'
+  []
+  [effective_strain_rate]
+    type = SR2Invariant
+    invariant_type = 'EFFECTIVE_STRAIN'
+    tensor = 'state/E_rate'
+    invariant = 'forces/effective_strain_rate'
+  []
+  [g]
+    type = KocksMeckingActivationEnergy
+    shear_modulus = 'mu'
+    k = 1.38064e-20
+    b = 2.474e-7
+    eps0 = 1e10
+  []
+  [flowrate]
+    type = KocksMeckingFlowSwitch
+    g0 = 'g0'
+    rate_independent_flow_rate = 'state/internal/gamma_rate_ri' 
+    rate_dependent_flow_rate = 'state/internal/gamma_rate_rd'
+    sharpness = 100.0
+  []
+  [Eprate]
+    type = AssociativePlasticFlow
+  []
+  [eprate]
+    type = AssociativeIsotropicPlasticHardening
+  []
+  [Erate]
+    type = SR2VariableRate
+    variable = 'state/E'
+    rate = 'state/E_rate'
+  []
+  [Eerate]
+    type = SR2LinearCombination
+    from_var = 'state/E_rate state/internal/Ep_rate'
+    to_var = 'state/internal/Ee_rate'
+    coefficients = '1 -1'
+  []
+  [elasticity]
+    type = LinearIsotropicElasticity
+    youngs_modulus = 1e5
+    poisson_ratio = 0.3
+    rate_form = true
+  []
+  [integrate_stress]
+    type = SR2BackwardEulerTimeIntegration
+    variable = 'state/S'
+  []
+  [integrate_ep]
+    type = ScalarBackwardEulerTimeIntegration
+    variable = 'state/internal/ep'
+  []
+  [mixed]
+    type = MixedControlSetup
+  []
+  [mixed_old]
+    type = MixedControlSetup
+    control = "old_forces/control"
+    mixed_state = "old_state/mixed_state"
+    fixed_values = "old_forces/fixed_values"
+    cauchy_stress = "old_state/S"
+    strain = "old_state/E"
+  []
+  [rename]
+    type = CopySR2
+    from = "residual/S"
+    to = "residual/mixed_state"
+  []
+  [implicit_rate]
+    type = ComposedModel
+    models = "isoharden elasticity g
+              mandel_stress vonmises
+              yield yield_zero normality eprate Eprate Erate Eerate
+              ri_flowrate rd_flowrate flowrate integrate_ep integrate_stress effective_strain_rate 
+              mixed mixed_old rename"
+  []
+[]