Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "LAROMANCEPartitionStressUpdateBase.h"
      11             : 
      12             : #include <libmesh/dense_matrix.h> // libMesh::cholesky_solve
      13             : 
      14             : registerMooseObjectAliased("TensorMechanicsApp",
      15             :                            LAROMANCEPartitionStressUpdateBase,
      16             :                            "LAROMANCEPartitionStressUpdate");
      17             : registerMooseObjectAliased("TensorMechanicsApp",
      18             :                            ADLAROMANCEPartitionStressUpdateBase,
      19             :                            "ADLAROMANCEPartitionStressUpdate");
      20             : 
      21             : template <bool is_ad>
      22             : InputParameters
      23          72 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::validParams()
      24             : {
      25          72 :   InputParameters params = LAROMANCEStressUpdateBaseTempl<is_ad>::validParams();
      26          72 :   params.addClassDescription("LAROMANCE base class for partitioned reduced order models");
      27         216 :   params.addRangeCheckedParam<Real>(
      28             :       "finite_difference_width",
      29         144 :       1.0e-2,
      30             :       "finite_difference_width > 0",
      31             :       "Factor multiplied against the input stress to compute the finite "
      32             :       "difference derivative of the parition weights.");
      33             : 
      34          72 :   return params;
      35           0 : }
      36             : 
      37             : template <bool is_ad>
      38          54 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::LAROMANCEPartitionStressUpdateBaseTempl(
      39             :     const InputParameters & parameters)
      40             :   : LAROMANCEStressUpdateBaseTempl<is_ad>(parameters),
      41         108 :     _finite_difference_width(this->template getParam<Real>("finite_difference_width"))
      42             : {
      43          54 : }
      44             : 
      45             : template <bool is_ad>
      46             : void
      47          54 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::initialSetup()
      48             : {
      49          54 :   LAROMANCEStressUpdateBaseTempl<is_ad>::initialSetup();
      50          54 :   _partition_Mmean = getClassificationMmean();
      51          54 :   _partition_Mscale = getClassificationMscale();
      52          54 :   _partition_Xu = getClassificationXu();
      53          54 :   _partition_Ell = getClassificationEll();
      54          54 :   _partition_Eta = getClassificationEta();
      55          54 :   _partition_Luu = getClassificationLuu();
      56          54 :   _partition_Vind = getClassificationVind();
      57          54 :   _partition_difference.resize(_partition_Xu[0].size(),
      58          54 :                                std::vector<GenericReal<is_ad>>(_partition_Xu.size()));
      59          54 :   _partition_distance.resize(_partition_difference.size());
      60          54 :   _partition_covariance.resize(_partition_distance.size());
      61          54 :   _partition_b.resize(_partition_covariance.size());
      62          54 :   _partition_A.resize(_partition_Luu[0].size(), _partition_Luu.size());
      63          54 : }
      64             : 
      65             : template <bool is_ad>
      66             : void
      67           9 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::exportJSON()
      68             : {
      69           9 :   LAROMANCEStressUpdateBaseTempl<is_ad>::exportJSON();
      70             : 
      71           9 :   this->_json["m_mean"] = getClassificationMmean();
      72           9 :   this->_json["m_scale"] = getClassificationMscale();
      73          18 :   this->_json["xu"] = getClassificationXu();
      74          18 :   this->_json["ell"] = getClassificationEll();
      75          18 :   this->_json["eta"] = getClassificationEta();
      76          18 :   this->_json["luu"] = getClassificationLuu();
      77           9 :   this->_json["vind"] = getClassificationVind().get_values();
      78           9 : }
      79             : 
      80             : template <bool is_ad>
      81             : void
      82        8232 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::computePartitionWeights(
      83             :     std::vector<GenericReal<is_ad>> & weights, std::vector<GenericReal<is_ad>> & dweights_dstress)
      84             : {
      85        8232 :   std::vector<bool> in_partition(_num_partitions, false);
      86        8232 :   unsigned int num_partitions_active = 0;
      87       24696 :   for (unsigned int p = 0; p < _num_partitions; ++p)
      88             :   {
      89       74088 :     for (unsigned int t = 0; t < this->_num_tiles[p]; ++t)
      90       57624 :       if (!in_partition[p] && this->checkInTile(p, t))
      91             :         in_partition[p] = true;
      92             : 
      93       16464 :     num_partitions_active += in_partition[p];
      94             :   }
      95             :   mooseAssert(num_partitions_active <= _num_partitions,
      96             :               "Number of paritions active must be less than total number of paritions");
      97        8232 :   if (!num_partitions_active)
      98           0 :     mooseException("Number of active partitions (",
      99             :                    num_partitions_active,
     100             :                    ") out of total number of partitions (",
     101             :                    _num_partitions,
     102             :                    ") is zero. This may be because you are trying to sample outside the total "
     103             :                    "applicability of the model");
     104             : 
     105        8232 :   if (_num_partitions == 1) // If only one parition present, all weights are one
     106             :   {
     107           0 :     weights[0] = 1.0;
     108           0 :     dweights_dstress[0] = 0.0;
     109             :   }
     110        8232 :   else if (_num_partitions == 2)
     111             :   {
     112        8232 :     if (num_partitions_active == 1) // If only present in one partition, math is easy
     113             :     {
     114           0 :       std::fill(dweights_dstress.begin(), dweights_dstress.end(), 0.0);
     115             : 
     116           0 :       if (in_partition[0])
     117           0 :         weights[1] = 0.0;
     118           0 :       else if (in_partition[1])
     119           0 :         weights[1] = 1.0;
     120             :       else
     121           0 :         mooseError("Internal error: Parition weight calculation incorrect when only one "
     122             :                    "partition is applicable");
     123             :     }
     124        8232 :     else if (num_partitions_active == 2) // If present in two partitions, compute weights
     125             :     {
     126        8232 :       weights[1] = computeSecondPartitionWeight();
     127        8232 :       computeDSecondPartitionWeightDStress(dweights_dstress[1]);
     128        8232 :       dweights_dstress[0] = -dweights_dstress[1];
     129             :     }
     130        8232 :     weights[0] = 1.0 - weights[1];
     131             :   }
     132             :   else
     133           0 :     mooseError("Internal error: number of partitions can only be 1 or 2");
     134             : 
     135        8232 :   if (_num_partitions == 2)
     136        8232 :     _second_partition_weight[_qp] = this->_partition_weights[1];
     137             :   else
     138           0 :     _second_partition_weight[_qp] = 0.0;
     139             : 
     140             :   // Given the weight, compute the sigmoid transformed value to smoothly transition
     141       24696 :   for (unsigned int p = 0; p < _num_partitions; ++p)
     142             :   {
     143       16464 :     weights[p] = (1.0 - sigmoid(0.0, 1.0, weights[p]));
     144       16464 :     dweights_dstress[p] *= -sigmoid(0.0, 1.0, weights[p], true);
     145             :   }
     146        8232 : }
     147             : 
     148             : template <bool is_ad>
     149             : GenericReal<is_ad>
     150       16464 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::computeSecondPartitionWeight()
     151             : {
     152             :   // extract and scale only the relevant inputs on which the GPR was trained (ignore
     153             :   // _old_strain_input_index)
     154       16464 :   std::vector<GenericReal<is_ad>> scaled_input_values(_num_inputs - 1);
     155             :   unsigned int inc = 0;
     156       98784 :   for (const auto i : index_range(_input_values))
     157             :   {
     158       82320 :     if (i != _old_strain_input_index)
     159             :     {
     160       65856 :       scaled_input_values[inc] =
     161       65856 :           (_input_values[i] - _partition_Mmean[inc]) / _partition_Mscale[inc];
     162       65856 :       inc++;
     163             :     }
     164             :   }
     165             : 
     166             :   // compute the distance of the new point to ALL training points and sum the distance over all
     167             :   // input variables first compute difference between points.
     168             :   // for-loop to compute difference
     169       82320 :   for (const auto i : index_range(_partition_Xu))
     170     3358656 :     for (const auto j : index_range(_partition_Xu[0]))
     171     3292800 :       _partition_difference[j][i] =
     172     3292800 :           Utility::pow<2>((_partition_Xu[i][j] - scaled_input_values[i]) / _partition_Ell);
     173             : 
     174             :   // sum difference over input dimension and take square root
     175           0 :   std::fill(_partition_distance.begin(), _partition_distance.end(), 0.0);
     176      839664 :   for (const auto i : index_range(_partition_difference))
     177             :   {
     178     4116000 :     for (const auto j : index_range(_partition_difference[0]))
     179     3292800 :       _partition_distance[i] += _partition_difference[i][j];
     180      823200 :     _partition_distance[i] = std::sqrt(_partition_distance[i]);
     181             :   }
     182             : 
     183             :   // compute the covariance wrt ALL training points: the larger the distance, the smaller the
     184             :   // covariance
     185      839664 :   for (const auto i : index_range(_partition_covariance))
     186      823200 :     _partition_covariance[i] =
     187      823200 :         Utility::pow<2>(_partition_Eta) * std::exp(-0.5 * _partition_distance[i]);
     188             : 
     189             :   // solve the system of equations: convert covariance to libMesh::DenseVector
     190      839664 :   for (const auto i : index_range(_partition_covariance))
     191      823200 :     _partition_b(i) = _partition_covariance[i];
     192             : 
     193             :   // convert Luu to libMesh::DenseMatrix, first fill values of A with Luu-values
     194      839664 :   for (const auto i : index_range(_partition_Luu[0]))
     195    41983200 :     for (const auto j : index_range(_partition_Luu))
     196    41160000 :       _partition_A(i, j) = _partition_Luu[j][i];
     197             : 
     198             :   // solve the linear system of equations (lu_solve)
     199       16464 :   DenseVector<GenericReal<is_ad>> ma(_partition_b.size());
     200       16464 :   _partition_A.lu_solve(_partition_b, ma);
     201             : 
     202             :   // induce full array through sparse grid and summarize into mu
     203       16464 :   GenericReal<is_ad> partition_weight = 0.0;
     204      839664 :   for (const auto i : index_range(ma))
     205      823200 :     partition_weight += ma(i) * _partition_Vind(i);
     206             : 
     207             :   // scale mu between 0 and 1: the weight of partition 2
     208       16464 :   partition_weight = -partition_weight + 1.0;
     209       16464 :   partition_weight = std::min(partition_weight, 1.0);
     210       16980 :   partition_weight = std::max(partition_weight, 0.0);
     211       16464 :   return partition_weight;
     212             : }
     213             : 
     214             : template <bool is_ad>
     215             : void
     216        8232 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::computeDSecondPartitionWeightDStress(
     217             :     GenericReal<is_ad> & dsecond_partition_weight_dstress)
     218             : {
     219        8232 :   const auto fd_tol =
     220        8232 :       _input_values[_stress_input_index] * _finite_difference_width; // Finite difference width
     221        8232 :   _input_values[_stress_input_index] += fd_tol;
     222        8232 :   dsecond_partition_weight_dstress =
     223        8232 :       (computeSecondPartitionWeight() - _second_partition_weight[_qp]) / fd_tol;
     224        8232 :   _input_values[_stress_input_index] -= fd_tol;
     225        8232 : }
     226             : 
     227             : template class LAROMANCEPartitionStressUpdateBaseTempl<false>;
     228             : template class LAROMANCEPartitionStressUpdateBaseTempl<true>;