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         144 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::validParams()
      24             : {
      25         144 :   InputParameters params = LAROMANCEStressUpdateBaseTempl<is_ad>::validParams();
      26         144 :   params.addClassDescription("LAROMANCE base class for partitioned reduced order models");
      27         432 :   params.addRangeCheckedParam<Real>(
      28             :       "finite_difference_width",
      29         288 :       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         144 :   return params;
      35           0 : }
      36             : 
      37             : template <bool is_ad>
      38         108 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::LAROMANCEPartitionStressUpdateBaseTempl(
      39             :     const InputParameters & parameters)
      40             :   : LAROMANCEStressUpdateBaseTempl<is_ad>(parameters),
      41         216 :     _finite_difference_width(this->template getParam<Real>("finite_difference_width"))
      42             : {
      43         108 : }
      44             : 
      45             : template <bool is_ad>
      46             : void
      47         108 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::initialSetup()
      48             : {
      49         108 :   LAROMANCEStressUpdateBaseTempl<is_ad>::initialSetup();
      50         108 :   _partition_Mmean = getClassificationMmean();
      51         108 :   _partition_Mscale = getClassificationMscale();
      52         108 :   _partition_Xu = getClassificationXu();
      53         108 :   _partition_Ell = getClassificationEll();
      54         108 :   _partition_Eta = getClassificationEta();
      55         108 :   _partition_Luu = getClassificationLuu();
      56         108 :   _partition_Vind = getClassificationVind();
      57         108 :   _partition_difference.resize(_partition_Xu[0].size(),
      58         108 :                                std::vector<GenericReal<is_ad>>(_partition_Xu.size()));
      59         108 :   _partition_distance.resize(_partition_difference.size());
      60         108 :   _partition_covariance.resize(_partition_distance.size());
      61         108 :   _partition_b.resize(_partition_covariance.size());
      62         108 :   _partition_A.resize(_partition_Luu[0].size(), _partition_Luu.size());
      63         108 : }
      64             : 
      65             : template <bool is_ad>
      66             : void
      67          18 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::exportJSON()
      68             : {
      69          18 :   LAROMANCEStressUpdateBaseTempl<is_ad>::exportJSON();
      70             : 
      71          18 :   this->_json["m_mean"] = getClassificationMmean();
      72          18 :   this->_json["m_scale"] = getClassificationMscale();
      73          36 :   this->_json["xu"] = getClassificationXu();
      74          36 :   this->_json["ell"] = getClassificationEll();
      75          36 :   this->_json["eta"] = getClassificationEta();
      76          36 :   this->_json["luu"] = getClassificationLuu();
      77          18 :   this->_json["vind"] = getClassificationVind().get_values();
      78          18 : }
      79             : 
      80             : template <bool is_ad>
      81             : void
      82       13616 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::computePartitionWeights(
      83             :     std::vector<GenericReal<is_ad>> & weights, std::vector<GenericReal<is_ad>> & dweights_dstress)
      84             : {
      85       13616 :   std::vector<bool> in_partition(_num_partitions, false);
      86       13616 :   unsigned int num_partitions_active = 0;
      87       40848 :   for (unsigned int p = 0; p < _num_partitions; ++p)
      88             :   {
      89      122544 :     for (unsigned int t = 0; t < this->_num_tiles[p]; ++t)
      90       95312 :       if (!in_partition[p] && this->checkInTile(p, t))
      91             :         in_partition[p] = true;
      92             : 
      93       27232 :     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       13616 :   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       13616 :   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       13616 :   else if (_num_partitions == 2)
     111             :   {
     112       13616 :     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       13616 :     else if (num_partitions_active == 2) // If present in two partitions, compute weights
     125             :     {
     126       13616 :       weights[1] = computeSecondPartitionWeight();
     127       13616 :       computeDSecondPartitionWeightDStress(dweights_dstress[1]);
     128       13616 :       dweights_dstress[0] = -dweights_dstress[1];
     129             :     }
     130       13616 :     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       13616 :   if (_num_partitions == 2)
     136       13616 :     _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       40848 :   for (unsigned int p = 0; p < _num_partitions; ++p)
     142             :   {
     143       27232 :     weights[p] = (1.0 - sigmoid(0.0, 1.0, weights[p]));
     144       27232 :     dweights_dstress[p] *= -sigmoid(0.0, 1.0, weights[p], true);
     145             :   }
     146       13616 : }
     147             : 
     148             : template <bool is_ad>
     149             : GenericReal<is_ad>
     150       27232 : 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       27232 :   std::vector<GenericReal<is_ad>> scaled_input_values(_num_inputs - 1);
     155             :   unsigned int inc = 0;
     156      163392 :   for (const auto i : index_range(_input_values))
     157             :   {
     158      136160 :     if (i != _old_strain_input_index)
     159             :     {
     160      108928 :       scaled_input_values[inc] =
     161      108928 :           (_input_values[i] - _partition_Mmean[inc]) / _partition_Mscale[inc];
     162      108928 :       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      136160 :   for (const auto i : index_range(_partition_Xu))
     170     5555328 :     for (const auto j : index_range(_partition_Xu[0]))
     171     5446400 :       _partition_difference[j][i] =
     172     5446400 :           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     1388832 :   for (const auto i : index_range(_partition_difference))
     177             :   {
     178     6808000 :     for (const auto j : index_range(_partition_difference[0]))
     179     5446400 :       _partition_distance[i] += _partition_difference[i][j];
     180     1361600 :     _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     1388832 :   for (const auto i : index_range(_partition_covariance))
     186     1361600 :     _partition_covariance[i] =
     187     1361600 :         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     1388832 :   for (const auto i : index_range(_partition_covariance))
     191     1361600 :     _partition_b(i) = _partition_covariance[i];
     192             : 
     193             :   // convert Luu to libMesh::DenseMatrix, first fill values of A with Luu-values
     194     1388832 :   for (const auto i : index_range(_partition_Luu[0]))
     195    69441600 :     for (const auto j : index_range(_partition_Luu))
     196    68080000 :       _partition_A(i, j) = _partition_Luu[j][i];
     197             : 
     198             :   // solve the linear system of equations (lu_solve)
     199       27232 :   DenseVector<GenericReal<is_ad>> ma(_partition_b.size());
     200       27232 :   _partition_A.lu_solve(_partition_b, ma);
     201             : 
     202             :   // induce full array through sparse grid and summarize into mu
     203       27232 :   GenericReal<is_ad> partition_weight = 0.0;
     204     1388832 :   for (const auto i : index_range(ma))
     205     1361600 :     partition_weight += ma(i) * _partition_Vind(i);
     206             : 
     207             :   // scale mu between 0 and 1: the weight of partition 2
     208       27232 :   partition_weight = -partition_weight + 1.0;
     209       27232 :   partition_weight = std::min(partition_weight, 1.0);
     210       28008 :   partition_weight = std::max(partition_weight, 0.0);
     211       27232 :   return partition_weight;
     212             : }
     213             : 
     214             : template <bool is_ad>
     215             : void
     216       13616 : LAROMANCEPartitionStressUpdateBaseTempl<is_ad>::computeDSecondPartitionWeightDStress(
     217             :     GenericReal<is_ad> & dsecond_partition_weight_dstress)
     218             : {
     219       13616 :   const auto fd_tol =
     220       13616 :       _input_values[_stress_input_index] * _finite_difference_width; // Finite difference width
     221       13616 :   _input_values[_stress_input_index] += fd_tol;
     222       13616 :   dsecond_partition_weight_dstress =
     223       13616 :       (computeSecondPartitionWeight() - _second_partition_weight[_qp]) / fd_tol;
     224       13616 :   _input_values[_stress_input_index] -= fd_tol;
     225       13616 : }
     226             : 
     227             : template class LAROMANCEPartitionStressUpdateBaseTempl<false>;
     228             : template class LAROMANCEPartitionStressUpdateBaseTempl<true>;