https://mooseframework.inl.gov
LAROMANCEPartitionStressUpdateBase.C
Go to the documentation of this file.
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 
11 
12 #include <libmesh/dense_matrix.h> // libMesh::cholesky_solve
13 
14 registerMooseObjectAliased("SolidMechanicsApp",
16  "LAROMANCEPartitionStressUpdate");
17 registerMooseObjectAliased("SolidMechanicsApp",
19  "ADLAROMANCEPartitionStressUpdate");
20 
21 template <bool is_ad>
24 {
26  params.addClassDescription("LAROMANCE base class for partitioned reduced order models");
27  params.addRangeCheckedParam<Real>(
28  "finite_difference_width",
29  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  return params;
35 }
36 
37 template <bool is_ad>
39  const InputParameters & parameters)
40  : LAROMANCEStressUpdateBaseTempl<is_ad>(parameters),
41  _finite_difference_width(this->template getParam<Real>("finite_difference_width"))
42 {
43 }
44 
45 template <bool is_ad>
46 void
48 {
50  _partition_Mmean = getClassificationMmean();
51  _partition_Mscale = getClassificationMscale();
52  _partition_Xu = getClassificationXu();
53  _partition_Ell = getClassificationEll();
54  _partition_Eta = getClassificationEta();
55  _partition_Luu = getClassificationLuu();
56  _partition_Vind = getClassificationVind();
57  _partition_difference.resize(_partition_Xu[0].size(),
58  std::vector<GenericReal<is_ad>>(_partition_Xu.size()));
59  _partition_distance.resize(_partition_difference.size());
60  _partition_covariance.resize(_partition_distance.size());
61  _partition_b.resize(_partition_covariance.size());
62  _partition_A.resize(_partition_Luu[0].size(), _partition_Luu.size());
63 }
64 
65 template <bool is_ad>
66 void
68 {
70 
71  this->_json["m_mean"] = getClassificationMmean();
72  this->_json["m_scale"] = getClassificationMscale();
73  this->_json["xu"] = getClassificationXu();
74  this->_json["ell"] = getClassificationEll();
75  this->_json["eta"] = getClassificationEta();
76  this->_json["luu"] = getClassificationLuu();
77  this->_json["vind"] = getClassificationVind().get_values();
78 }
79 
80 template <bool is_ad>
81 void
83  std::vector<GenericReal<is_ad>> & weights, std::vector<GenericReal<is_ad>> & dweights_dstress)
84 {
85  std::vector<bool> in_partition(_num_partitions, false);
86  unsigned int num_partitions_active = 0;
87  for (unsigned int p = 0; p < _num_partitions; ++p)
88  {
89  for (unsigned int t = 0; t < this->_num_tiles[p]; ++t)
90  if (!in_partition[p] && this->checkInTile(p, t))
91  in_partition[p] = true;
92 
93  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  if (!num_partitions_active)
98  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  if (_num_partitions == 1) // If only one parition present, all weights are one
106  {
107  weights[0] = 1.0;
108  dweights_dstress[0] = 0.0;
109  }
110  else if (_num_partitions == 2)
111  {
112  if (num_partitions_active == 1) // If only present in one partition, math is easy
113  {
114  std::fill(dweights_dstress.begin(), dweights_dstress.end(), 0.0);
115 
116  if (in_partition[0])
117  weights[1] = 0.0;
118  else if (in_partition[1])
119  weights[1] = 1.0;
120  else
121  mooseError("Internal error: Parition weight calculation incorrect when only one "
122  "partition is applicable");
123  }
124  else if (num_partitions_active == 2) // If present in two partitions, compute weights
125  {
126  weights[1] = computeSecondPartitionWeight();
127  computeDSecondPartitionWeightDStress(dweights_dstress[1]);
128  dweights_dstress[0] = -dweights_dstress[1];
129  }
130  weights[0] = 1.0 - weights[1];
131  }
132  else
133  mooseError("Internal error: number of partitions can only be 1 or 2");
134 
135  if (_num_partitions == 2)
136  _second_partition_weight[_qp] = this->_partition_weights[1];
137  else
138  _second_partition_weight[_qp] = 0.0;
139 
140  // Given the weight, compute the sigmoid transformed value to smoothly transition
141  for (unsigned int p = 0; p < _num_partitions; ++p)
142  {
143  weights[p] = (1.0 - sigmoid(0.0, 1.0, weights[p]));
144  dweights_dstress[p] *= -sigmoid(0.0, 1.0, weights[p], true);
145  }
146 }
147 
148 template <bool is_ad>
151 {
152  // extract and scale only the relevant inputs on which the GPR was trained (ignore
153  // _old_strain_input_index)
154  std::vector<GenericReal<is_ad>> scaled_input_values(_num_inputs - 1);
155  unsigned int inc = 0;
156  for (const auto i : index_range(_input_values))
157  {
158  if (i != _old_strain_input_index)
159  {
160  scaled_input_values[inc] =
161  (_input_values[i] - _partition_Mmean[inc]) / _partition_Mscale[inc];
162  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  for (const auto i : index_range(_partition_Xu))
170  for (const auto j : index_range(_partition_Xu[0]))
171  _partition_difference[j][i] =
172  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  std::fill(_partition_distance.begin(), _partition_distance.end(), 0.0);
176  for (const auto i : index_range(_partition_difference))
177  {
178  for (const auto j : index_range(_partition_difference[0]))
179  _partition_distance[i] += _partition_difference[i][j];
180  _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  for (const auto i : index_range(_partition_covariance))
186  _partition_covariance[i] =
187  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  for (const auto i : index_range(_partition_covariance))
191  _partition_b(i) = _partition_covariance[i];
192 
193  // convert Luu to libMesh::DenseMatrix, first fill values of A with Luu-values
194  for (const auto i : index_range(_partition_Luu[0]))
195  for (const auto j : index_range(_partition_Luu))
196  _partition_A(i, j) = _partition_Luu[j][i];
197 
198  // solve the linear system of equations (lu_solve)
199  DenseVector<GenericReal<is_ad>> ma(_partition_b.size());
200  _partition_A.lu_solve(_partition_b, ma);
201 
202  // induce full array through sparse grid and summarize into mu
203  GenericReal<is_ad> partition_weight = 0.0;
204  for (const auto i : index_range(ma))
205  partition_weight += ma(i) * _partition_Vind(i);
206 
207  // scale mu between 0 and 1: the weight of partition 2
208  partition_weight = -partition_weight + 1.0;
209  partition_weight = std::min(partition_weight, 1.0);
210  partition_weight = std::max(partition_weight, 0.0);
211  return partition_weight;
212 }
213 
214 template <bool is_ad>
215 void
217  GenericReal<is_ad> & dsecond_partition_weight_dstress)
218 {
219  const auto fd_tol =
220  _input_values[_stress_input_index] * _finite_difference_width; // Finite difference width
221  _input_values[_stress_input_index] += fd_tol;
222  dsecond_partition_weight_dstress =
223  (computeSecondPartitionWeight() - _second_partition_weight[_qp]) / fd_tol;
224  _input_values[_stress_input_index] -= fd_tol;
225 }
226 
Moose::GenericType< Real, is_ad > GenericReal
virtual void computePartitionWeights(std::vector< GenericReal< is_ad >> &weights, std::vector< GenericReal< is_ad >> &dweights_dstress) override
Compute the weight of the different partitions.
void mooseError(Args &&... args)
virtual GenericReal< is_ad > computeSecondPartitionWeight()
Compute the partition weight on the location in input-space, based on a calibrated Gaussian Process R...
virtual void computeDSecondPartitionWeightDStress(GenericReal< is_ad > &dsecond_partition_weight_dstress)
Compute the derivative of the partition weight of the second partition w.r.t.
LAROMANCEPartitionStressUpdateBaseTempl(const InputParameters &parameters)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void addClassDescription(const std::string &doc_string)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
auto index_range(const T &sizable)
registerMooseObjectAliased("SolidMechanicsApp", LAROMANCEPartitionStressUpdateBase, "LAROMANCEPartitionStressUpdate")