Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       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("SolidMechanicsApp",
      15             :                            LAROMANCEPartitionStressUpdateBase,
      16             :                            "LAROMANCEPartitionStressUpdate");
      17             : registerMooseObjectAliased("SolidMechanicsApp",
      18             :                            ADLAROMANCEPartitionStressUpdateBase,
      19             :                            "ADLAROMANCEPartitionStressUpdate");
      20             : 
      21             : template <bool is_ad>
      22             : InputParameters
      23          96 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::validParams()
      24             : {
      25          96 :   InputParameters params = LAROMANCEStressUpdateBaseTempl<is_ad>::validParams();
      26          96 :   params.addClassDescription("LAROMANCE base class for partitioned reduced order models");
      27         288 :   params.addRangeCheckedParam<Real>(
      28             :       "finite_difference_width",
      29         192 :       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          96 :   return params;
      35           0 : }
      36             : 
      37             : template <bool is_ad>
      38          72 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::LAROMANCEPartitionStressUpdateBaseTempl(
      39             :     const InputParameters & parameters)
      40             :   : LAROMANCEStressUpdateBaseTempl<is_ad>(parameters),
      41         144 :     _finite_difference_width(this->template getParam<Real>("finite_difference_width"))
      42             : {
      43          72 : }
      44             : 
      45             : template <bool is_ad>
      46             : void
      47          72 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::initialSetup()
      48             : {
      49          72 :   LAROMANCEStressUpdateBaseTempl<is_ad>::initialSetup();
      50          72 :   _partition_Mmean = getClassificationMmean();
      51          72 :   _partition_Mscale = getClassificationMscale();
      52          72 :   _partition_Xu = getClassificationXu();
      53          72 :   _partition_Ell = getClassificationEll();
      54          72 :   _partition_Eta = getClassificationEta();
      55          72 :   _partition_Luu = getClassificationLuu();
      56          72 :   _partition_Vind = getClassificationVind();
      57          72 :   _partition_difference.resize(_partition_Xu[0].size(),
      58         144 :                                std::vector<GenericReal<is_ad>>(_partition_Xu.size()));
      59          72 :   _partition_distance.resize(_partition_difference.size());
      60          72 :   _partition_covariance.resize(_partition_distance.size());
      61          72 :   _partition_b.resize(_partition_covariance.size());
      62          72 :   _partition_A.resize(_partition_Luu[0].size(), _partition_Luu.size());
      63          72 : }
      64             : 
      65             : template <bool is_ad>
      66             : void
      67          12 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::exportJSON()
      68             : {
      69          12 :   LAROMANCEStressUpdateBaseTempl<is_ad>::exportJSON();
      70             : 
      71          24 :   this->_json["m_mean"] = getClassificationMmean();
      72          24 :   this->_json["m_scale"] = getClassificationMscale();
      73          24 :   this->_json["xu"] = getClassificationXu();
      74          24 :   this->_json["ell"] = getClassificationEll();
      75          24 :   this->_json["eta"] = getClassificationEta();
      76          24 :   this->_json["luu"] = getClassificationLuu();
      77          12 :   this->_json["vind"] = getClassificationVind().get_values();
      78          12 : }
      79             : 
      80             : template <bool is_ad>
      81             : void
      82       10813 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::computePartitionWeights(
      83             :     std::vector<GenericReal<is_ad>> & weights, std::vector<GenericReal<is_ad>> & dweights_dstress)
      84             : {
      85       10813 :   std::vector<bool> in_partition(_num_partitions, false);
      86       10813 :   unsigned int num_partitions_active = 0;
      87       32439 :   for (unsigned int p = 0; p < _num_partitions; ++p)
      88             :   {
      89       97317 :     for (unsigned int t = 0; t < this->_num_tiles[p]; ++t)
      90       75691 :       if (!in_partition[p] && this->checkInTile(p, t))
      91             :         in_partition[p] = true;
      92             : 
      93       21626 :     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       10813 :   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       10813 :   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       10813 :   else if (_num_partitions == 2)
     111             :   {
     112       10813 :     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       10813 :     else if (num_partitions_active == 2) // If present in two partitions, compute weights
     125             :     {
     126       10813 :       weights[1] = computeSecondPartitionWeight();
     127       10813 :       computeDSecondPartitionWeightDStress(dweights_dstress[1]);
     128       10813 :       dweights_dstress[0] = -dweights_dstress[1];
     129             :     }
     130       10813 :     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       10813 :   if (_num_partitions == 2)
     136       10813 :     _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       32439 :   for (unsigned int p = 0; p < _num_partitions; ++p)
     142             :   {
     143       21626 :     weights[p] = (1.0 - sigmoid(0.0, 1.0, weights[p]));
     144       21626 :     dweights_dstress[p] *= -sigmoid(0.0, 1.0, weights[p], true);
     145             :   }
     146       10813 : }
     147             : 
     148             : template <bool is_ad>
     149             : GenericReal<is_ad>
     150       21626 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::computeSecondPartitionWeight()
     151             : {
     152             :   using std::sqrt, std::exp, std::min, std::max;
     153             : 
     154             :   // extract and scale only the relevant inputs on which the GPR was trained (ignore
     155             :   // _old_strain_input_index)
     156       21626 :   std::vector<GenericReal<is_ad>> scaled_input_values(_num_inputs - 1);
     157             :   unsigned int inc = 0;
     158      129756 :   for (const auto i : index_range(_input_values))
     159             :   {
     160      108130 :     if (i != _old_strain_input_index)
     161             :     {
     162       86504 :       scaled_input_values[inc] =
     163       86504 :           (_input_values[i] - _partition_Mmean[inc]) / _partition_Mscale[inc];
     164       86504 :       inc++;
     165             :     }
     166             :   }
     167             : 
     168             :   // compute the distance of the new point to ALL training points and sum the distance over all
     169             :   // input variables first compute difference between points.
     170             :   // for-loop to compute difference
     171      108130 :   for (const auto i : index_range(_partition_Xu))
     172     4411704 :     for (const auto j : index_range(_partition_Xu[0]))
     173     4325200 :       _partition_difference[j][i] =
     174     4325200 :           Utility::pow<2>((_partition_Xu[i][j] - scaled_input_values[i]) / _partition_Ell);
     175             : 
     176             :   // sum difference over input dimension and take square root
     177           0 :   std::fill(_partition_distance.begin(), _partition_distance.end(), 0.0);
     178     1102926 :   for (const auto i : index_range(_partition_difference))
     179             :   {
     180     5406500 :     for (const auto j : index_range(_partition_difference[0]))
     181     4325200 :       _partition_distance[i] += _partition_difference[i][j];
     182     1081300 :     _partition_distance[i] = sqrt(_partition_distance[i]);
     183             :   }
     184             : 
     185             :   // compute the covariance wrt ALL training points: the larger the distance, the smaller the
     186             :   // covariance
     187     1102926 :   for (const auto i : index_range(_partition_covariance))
     188     1081300 :     _partition_covariance[i] = Utility::pow<2>(_partition_Eta) * exp(-0.5 * _partition_distance[i]);
     189             : 
     190             :   // solve the system of equations: convert covariance to libMesh::DenseVector
     191     1102926 :   for (const auto i : index_range(_partition_covariance))
     192     1081300 :     _partition_b(i) = _partition_covariance[i];
     193             : 
     194             :   // convert Luu to libMesh::DenseMatrix, first fill values of A with Luu-values
     195     1102926 :   for (const auto i : index_range(_partition_Luu[0]))
     196    55146300 :     for (const auto j : index_range(_partition_Luu))
     197    54065000 :       _partition_A(i, j) = _partition_Luu[j][i];
     198             : 
     199             :   // solve the linear system of equations (lu_solve)
     200       21626 :   DenseVector<GenericReal<is_ad>> ma(_partition_b.size());
     201       21626 :   _partition_A.lu_solve(_partition_b, ma);
     202             : 
     203             :   // induce full array through sparse grid and summarize into mu
     204       21626 :   GenericReal<is_ad> partition_weight = 0.0;
     205     1102926 :   for (const auto i : index_range(ma))
     206     1081300 :     partition_weight += ma(i) * _partition_Vind(i);
     207             : 
     208             :   // scale mu between 0 and 1: the weight of partition 2
     209       21626 :   partition_weight = -partition_weight + 1.0;
     210       21626 :   partition_weight = min(partition_weight, 1.0);
     211       22244 :   partition_weight = max(partition_weight, 0.0);
     212       21626 :   return partition_weight;
     213       21626 : }
     214             : 
     215             : template <bool is_ad>
     216             : void
     217       10813 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::computeDSecondPartitionWeightDStress(
     218             :     GenericReal<is_ad> & dsecond_partition_weight_dstress)
     219             : {
     220       10813 :   const auto fd_tol =
     221       10813 :       _input_values[_stress_input_index] * _finite_difference_width; // Finite difference width
     222       10813 :   _input_values[_stress_input_index] += fd_tol;
     223       10813 :   dsecond_partition_weight_dstress =
     224       10813 :       (computeSecondPartitionWeight() - _second_partition_weight[_qp]) / fd_tol;
     225       10813 :   _input_values[_stress_input_index] -= fd_tol;
     226       10813 : }
     227             : 
     228             : template class LAROMANCEPartitionStressUpdateBaseTempl<false>;
     229             : template class LAROMANCEPartitionStressUpdateBaseTempl<true>;