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 "ParameterMeshFunction.h"
      11             : 
      12             : #include "AddVariableAction.h"
      13             : 
      14             : registerMooseObject("OptimizationApp", ParameterMeshFunction);
      15             : 
      16             : InputParameters
      17         636 : ParameterMeshFunction::validParams()
      18             : {
      19         636 :   InputParameters params = OptimizationFunction::validParams();
      20         636 :   params.addClassDescription("Optimization function with parameters represented by a mesh and "
      21             :                              "finite-element shape functions.");
      22             : 
      23        1272 :   params.addRequiredParam<FileName>("exodus_mesh", "File containing parameter mesh.");
      24        1272 :   params.addParam<MooseEnum>("family",
      25        1272 :                              AddVariableAction::getNonlinearVariableFamilies(),
      26             :                              "Family of FE shape functions for parameter.");
      27        1272 :   params.addParam<MooseEnum>("order",
      28        1272 :                              AddVariableAction::getNonlinearVariableOrders(),
      29             :                              "Order of FE shape functions for parameter.");
      30             : 
      31        1272 :   params.addRequiredParam<ReporterName>(
      32             :       "parameter_name", "Reporter or VectorPostprocessor vector containing parameter values.");
      33        1272 :   params.addParam<ReporterName>("time_name",
      34             :                                 "Name of vector-postprocessor or reporter vector containing time, "
      35             :                                 "default assumes time independence.");
      36        1272 :   params.addParam<bool>(
      37        1272 :       "project_points", false, "Whether to find the closest point on parameter mesh.");
      38        1272 :   params.addParam<unsigned int>(
      39             :       "kdtree_candidates",
      40        1272 :       5,
      41             :       "Number of nearest node candidates to consider when projecting points to parameter mesh.");
      42             : 
      43         636 :   return params;
      44           0 : }
      45             : 
      46         334 : ParameterMeshFunction::ParameterMeshFunction(const InputParameters & parameters)
      47             :   : OptimizationFunction(parameters),
      48             :     ReporterInterface(this),
      49         666 :     _parameter_mesh(AddVariableAction::feType(parameters),
      50         666 :                     getParam<FileName>("exodus_mesh"),
      51             :                     {},
      52        1000 :                     getParam<bool>("project_points"),
      53         666 :                     getParam<unsigned int>("kdtree_candidates")),
      54         332 :     _values(getReporterValue<std::vector<Real>>("parameter_name")),
      55         686 :     _coordt(isParamValid("time_name") ? getReporterValue<std::vector<Real>>("time_name")
      56         334 :                                       : _empty_vec)
      57             : {
      58         332 : }
      59             : 
      60             : Real
      61     6337717 : ParameterMeshFunction::value(Real t, const Point & p) const
      62             : {
      63     6337717 :   checkSize();
      64             : 
      65     6337715 :   const auto ti = interpolateTime(t);
      66     6337715 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
      67     6337715 :   const dof_id_type offset1 = ti[1].first * _parameter_mesh.size();
      68             : 
      69             :   std::vector<dof_id_type> dof_indices;
      70             :   std::vector<Real> weights;
      71             : 
      72     6337715 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
      73             : 
      74             :   Real val = 0;
      75    32223880 :   for (const auto & i : index_range(dof_indices))
      76    25886165 :     val += (_values[dof_indices[i] + offset0] * ti[0].second +
      77    25886165 :             _values[dof_indices[i] + offset1] * ti[1].second) *
      78             :            weights[i];
      79     6337715 :   return val;
      80     6337715 : }
      81             : 
      82             : RealGradient
      83       16200 : ParameterMeshFunction::gradient(Real t, const Point & p) const
      84             : {
      85       16200 :   checkSize();
      86             : 
      87       16200 :   const auto ti = interpolateTime(t);
      88       16200 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
      89       16200 :   const dof_id_type offset1 = ti[1].first * _parameter_mesh.size();
      90             : 
      91             :   std::vector<dof_id_type> dof_indices;
      92             :   std::vector<RealGradient> weights;
      93       16200 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
      94             : 
      95             :   RealGradient val(0, 0, 0);
      96       84600 :   for (const auto & i : index_range(dof_indices))
      97       68400 :     val += (_values[dof_indices[i] + offset0] * ti[0].second +
      98       68400 :             _values[dof_indices[i] + offset1] * ti[1].second) *
      99             :            weights[i];
     100       16200 :   return val;
     101       16200 : }
     102             : 
     103             : Real
     104           0 : ParameterMeshFunction::timeDerivative(Real t, const Point & p) const
     105             : {
     106           0 :   checkSize();
     107             : 
     108           0 :   const auto ti = interpolateTime(t, true);
     109           0 :   if (ti[0].first == ti[1].first)
     110             :     return 0.0;
     111           0 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
     112           0 :   const dof_id_type offset1 = ti[1].first * _parameter_mesh.size();
     113             : 
     114             :   std::vector<dof_id_type> dof_indices;
     115             :   std::vector<Real> weights;
     116           0 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
     117             : 
     118             :   Real val = 0;
     119           0 :   for (const auto & i : index_range(dof_indices))
     120           0 :     val += (_values[dof_indices[i] + offset0] * ti[0].second +
     121           0 :             _values[dof_indices[i] + offset1] * ti[1].second) *
     122             :            weights[i];
     123             :   return val;
     124           0 : }
     125             : 
     126             : std::vector<Real>
     127      598601 : ParameterMeshFunction::parameterGradient(Real t, const Point & p) const
     128             : {
     129      598601 :   const auto ti = interpolateTime(t);
     130      598601 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
     131      598601 :   const dof_id_type offset1 = ti[1].first * _parameter_mesh.size();
     132             : 
     133             :   std::vector<dof_id_type> dof_indices;
     134             :   std::vector<Real> weights;
     135      598601 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
     136             : 
     137             :   dof_id_type sz = _parameter_mesh.size();
     138      598601 :   if (!_coordt.empty())
     139        6534 :     sz *= _coordt.size();
     140      598601 :   std::vector<Real> pg(sz, 0.0);
     141     3054383 :   for (const auto & i : index_range(dof_indices))
     142             :   {
     143     2455782 :     pg[dof_indices[i] + offset0] += weights[i] * ti[0].second;
     144     2455782 :     pg[dof_indices[i] + offset1] += weights[i] * ti[1].second;
     145             :   }
     146      598601 :   return pg;
     147      598601 : }
     148             : 
     149             : std::array<std::pair<std::size_t, Real>, 2>
     150     6952516 : ParameterMeshFunction::interpolateTime(Real t, bool derivative) const
     151             : {
     152     6952516 :   std::array<std::pair<std::size_t, Real>, 2> ti;
     153     6952516 :   if (_coordt.size() <= 1 || MooseUtils::absoluteFuzzyLessEqual(t, _coordt[0]))
     154             :   {
     155             :     ti[0] = {0, 0.5};
     156             :     ti[1] = {0, 0.5};
     157             :   }
     158       18945 :   else if (MooseUtils::absoluteFuzzyGreaterEqual(t, _coordt.back()))
     159             :   {
     160        3789 :     ti[0] = {_coordt.size() - 1, 0.5};
     161             :     ti[1] = {_coordt.size() - 1, 0.5};
     162             :   }
     163             :   else
     164       22734 :     for (std::size_t i = 1; i < _coordt.size(); ++i)
     165       22734 :       if (MooseUtils::absoluteFuzzyGreaterEqual(_coordt[i], t))
     166             :       {
     167       15156 :         const Real dt = _coordt[i] - _coordt[i - 1];
     168       15156 :         if (derivative)
     169             :         {
     170           0 :           ti[0] = {i - 1, -1.0 / dt};
     171             :           ti[1] = {i, 1.0 / dt};
     172             :         }
     173             :         else
     174             :         {
     175       15156 :           ti[0] = {i - 1, (_coordt[i] - t) / dt};
     176       15156 :           ti[1] = {i, (t - _coordt[i - 1]) / dt};
     177             :         }
     178             :         break;
     179             :       }
     180             : 
     181     6952516 :   return ti;
     182             : }
     183             : 
     184             : void
     185     6353917 : ParameterMeshFunction::checkSize() const
     186             : {
     187             :   dof_id_type sz = _parameter_mesh.size();
     188     6353917 :   if (!_coordt.empty())
     189       16200 :     sz *= _coordt.size();
     190     6353917 :   if (sz != _values.size())
     191           2 :     paramError("parameter_name",
     192             :                "Size of parameter vector (",
     193             :                _values.size(),
     194             :                ") does not match number of degrees of freedom in mesh (",
     195             :                sz,
     196             :                ").");
     197     6353915 : }